# Copyright (C) 2022 VRVis.
# All rights reserved.
# Contact: VRVis Forschungs-GmbH (office@vrvis.at)
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# 1. Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
# 2. Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
# 3. All advertising materials mentioning features or use of this software
# must display the following acknowledgement:
# This product includes software developed by the VRVis Forschungs-GmbH.
# 4. Neither the name of the VRVis Forschungs-GmbH nor the
# names of its contributors may be used to endorse or promote products
# derived from this software without specific prior written permission.
#
#setup
workingDir<-"C:/Users/ganglberger/Documents/Evolution/evo_paper_v3_olga_data/R_MAIN_CODE"
knitr::opts_knit$set(root.dir = normalizePath(workingDir))
setwd(workingDir)
set.seed(42)
#load libraries
#Runs under R 4.0.3
options(java.parameters = "- Xmx2024m")
library("xlsx")
library("corrplot")
## corrplot 0.84 loaded
library("R.matlab")
## R.matlab v3.6.2 (2018-09-26) successfully loaded. See ?R.matlab for help.
##
## Attaching package: 'R.matlab'
## The following objects are masked from 'package:base':
##
## getOption, isOpen
library("rjson")
library("gplots")
##
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
##
## lowess
library("png")
library("e1071")
library("rsvg")
library("XML")
library("igraph")
##
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
## The following object is masked from 'package:base':
##
## union
library('knitr')
library("biomaRt")
library('org.Hs.eg.db')
## Loading required package: AnnotationDbi
## Loading required package: stats4
## Loading required package: BiocGenerics
## Loading required package: parallel
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:igraph':
##
## normalize, path, union
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which, which.max, which.min
## Loading required package: Biobase
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: IRanges
## Loading required package: S4Vectors
##
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:gplots':
##
## space
## The following object is masked from 'package:base':
##
## expand.grid
##
## Attaching package: 'IRanges'
## The following object is masked from 'package:grDevices':
##
## windows
##
#Read DNDS-ratio data
dnds<-read.csv2(dataFile,stringsAsFactors=FALSE,dec = ".",sep=",",header=length(colsToUseNames)==0)
# values below zero mean:
#-1 - unstable dn/ds in 5 runs
#-2 - short edge, unstable
#-3 - 0 values, dN=0
#-4 - 999 values, dS=0
#-5 identical sequences (the majority between HUMAN and NEANDERTHAL)
if(length(colsToUseNames)>0){
colnames(dnds)[colsToUse]<-colsToUseNames
}
humanEntrez<-dnds[,2]
print(paste0("Dimension of data: ",dim(dnds)[1]," x ",dim(dnds)[2]))
## [1] "Dimension of data: 8978 x 25"
head(dnds)
#Read ontology (hierarchy of regions + a color palette for the output of MRI style plots)
documentHuman <- fromJSON(file="../storage/ontologyHuman.json", method='C')
#gets acronym of a region given by a region id (from ontology file)
getAcronymByID <- function(child,ID){
if(!is.null(child)){
if(child$id==ID){
return(child)
}
}
for(actchild in child$children){
res<-getAcronymByID(actchild,ID)
if(!is.null(res)){
if(res$id==ID){
return(res)
}
}
}
return(NULL)
}
getIDbyName <- function(child,name){
if(!is.null(child)){
if(child$name==name){
return(child)
}
}
for(actchild in child$children){
res<-getIDbyName(actchild,name)
if(!is.null(res)){
if(res$name==name){
return(res)
}
}
}
return(NULL)
}
jet.colors <-
colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
"#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
#Read Human data and filter it for genes that are also in the DNDS-ratio file
#gives a vector with the size of the atlas (1D vector of human brain atlas biopsy sites) where every position
#is zero, execpt for positions where the region (given by the ID) can be found
getAtlasRegionsOfIDHuman <- function(atlasRegions,ID){
newAtlasRegions<-atlasRegions==ID
childrens<-getAcronymByID(documentHuman$msg[[1]],ID)$children
while(length(childrens)>0){
newChildren<-c()
for(i in 1:length(childrens)){
newAtlasRegions<-newAtlasRegions|(atlasRegions==childrens[[i]]$id)
newChildren<-c(newChildren,getAcronymByID(documentHuman$msg[[1]],childrens[[i]]$id)$children)
}
childrens<-newChildren
}
return(newAtlasRegions)
}
atlasRegionsHuman<-readMat("../storage/atlasRegionsHuman.mat") #1D vector of human brain atlas biopsy sites
regionInfo<-cbind(atlasRegionsHuman$indexX,atlasRegionsHuman$indexY,atlasRegionsHuman$indexZ)
atlasRegionsHuman<-as.vector(atlasRegionsHuman$atlasRegions)
regionName<-c()
regionShortName<-c()
for(rid in atlasRegionsHuman){
a<-getAcronymByID(documentHuman$msg[[1]],rid)
regionName<-c(regionName,a$name)
regionShortName<-c(regionShortName,a$acronym)
}
regionInfo<-cbind(regionShortName,regionName,regionInfo)
colnames(regionInfo)<-c("acronym","name","mniX","mniY","mniZ")
if(!exists("expressionHuman") && !exists("entrezIDsOfRowsHuman")){ #this if just makes it faster to reload biopsy site level gene expression data, since readMat is quite slow
all_genes<-readMat("../storage/all_genes_expression_Human.mat")
expressionHuman<-all_genes$expressionMatrix
entrezIDsOfRowsHuman<-all_genes$entrezIDsOfRows
remove(all_genes)
gc()
expressionHuman<-t(apply(expressionHuman,1,function(x){(x-median(x,na.rm=TRUE))/mad(x,na.rm=TRUE)})) #normalize robust per biopsy site
expressionHuman<-apply(expressionHuman,2,function(x){(x-median(x,na.rm=TRUE))/mad(x,na.rm=TRUE)}) #normalize robust per gene (i.e. over the brain)
}
##
## Attaching package: 'Matrix'
## The following object is masked from 'package:S4Vectors':
##
## expand
print(paste0("Dimension of human data: ",dim(expressionHuman)[1]," x ",dim(expressionHuman)[2]))
## [1] "Dimension of human data: 3702 x 20787"
entrezIDsInAMBAandDNDSHuman<-intersect(unique(humanEntrez),unique(entrezIDsOfRowsHuman))
print(paste0("Human genes that are in ABA and DNDS-ratio file: ",length(entrezIDsInAMBAandDNDSHuman)))
## [1] "Human genes that are in ABA and DNDS-ratio file: 8977"
#Filter genes that are human (and are in the DNDS ratio file)
humanEntrezInAmba<-unique(humanEntrez[is.element(humanEntrez,entrezIDsInAMBAandDNDSHuman)])
print(paste0("Genes that are in human (and are in the DNDS ratio file): ",length(humanEntrezInAmba), " unique: ",length(unique(humanEntrezInAmba))))
## [1] "Genes that are in human (and are in the DNDS ratio file): 8977 unique: 8977"
#Select columns (branches) that should be used. Take only rows (genes) that have no NAs in the branches and are not all 0
#Take also the mean for genes that are not unique for a human
dnds[,colsToUse]<-apply(dnds[,colsToUse],2,function(x){as.numeric(x)})
dnds_org<-dnds
dnds[,colsToUse][(dnds[,colsToUse]==(-3)) & !is.na(dnds[,colsToUse])]<-0 #set DNDS ratios wih a -3 value (dN=0, dS>0) to 0
dnds[,colsToUse][dnds[,colsToUse]<0 & !is.na(dnds[,colsToUse])]<-NA # set other special values below 0 to NA
dnds_org[,colsToUse][dnds_org[,colsToUse]<0 & !is.na(dnds_org[,colsToUse])]<-NA
print(paste0("Use cols: "))
## [1] "Use cols: "
print(paste0(names(dnds)[colsToUse]))
## [1] "MOUSE" "OTOGA" "CALJA" "MACMU" "NOMLE"
## [6] "GORGO" "PANTR" "DENISOVAN" "NEANDERTHAL" "HUMAN"
## [11] "MO_CMNGPDNH" "MOC_MNGPDNH" "MOCM_NGPDNH" "MOCMN_GPDNH" "MOCMNG_HNDP"
## [16] "MOCMNGP_HND" "MOCMNGPD_HN" "MOCMNP_HNDG" "MOCMNDHN_PG" "MOCMNGPH_ND"
## [21] "MOCMNGPN_HD"
print(paste0("Rows before filter: ",dim(dnds)[1]))
## [1] "Rows before filter: 8978"
dnds_per_gene<-matrix(0,nrow=length(humanEntrezInAmba),ncol=length(colsToUse))
for(actCol in 1:length(colsToUse)){
for(actRow in 1:length(humanEntrezInAmba)){
dnds_per_gene[actRow,actCol]<-mean(dnds[humanEntrez==humanEntrezInAmba[actRow],colsToUse[actCol]],na.rm=TRUE)
}
}
dnds_per_gene[is.nan(dnds_per_gene)]<-NA
print(paste0("Rows after filter for unique human of the selected columns: ",dim(dnds_per_gene)[1]))
## [1] "Rows after filter for unique human of the selected columns: 8977"
rownames(dnds_per_gene)<-humanEntrezInAmba
colnames(dnds_per_gene)<-names(dnds)[colsToUse]
print(paste0("Dimension of data: ",dim(dnds_per_gene)[1]," x ",dim(dnds_per_gene)[2]," length of entrez (should be first dimension: ",length(humanEntrezInAmba)))
## [1] "Dimension of data: 8977 x 21 length of entrez (should be first dimension: 8977"
print("Amount of rows not NA")
## [1] "Amount of rows not NA"
apply(as.matrix(dnds_per_gene),2,function(x){sum(!is.na(x))})
## MOUSE OTOGA CALJA MACMU NOMLE GORGO
## 8936 8917 8865 8832 8737 7875
## PANTR DENISOVAN NEANDERTHAL HUMAN MO_CMNGPDNH MOC_MNGPDNH
## 7261 8836 1097 2222 8786 8350
## MOCM_NGPDNH MOCMN_GPDNH MOCMNG_HNDP MOCMNGP_HND MOCMNGPD_HN MOCMNP_HNDG
## 8013 8149 3899 6893 1199 639
## MOCMNDHN_PG MOCMNGPH_ND MOCMNGPN_HD
## 652 1507 54
#Rank-Normalize dnds column wise
print(paste0("Column means/sd before normalization: "))
## [1] "Column means/sd before normalization: "
cbind(apply(as.matrix(dnds_per_gene),2,function(x){mean(x,na.rm=TRUE)}),apply(as.matrix(dnds_per_gene),2,function(x){sd(x,na.rm=TRUE)}))
## [,1] [,2]
## MOUSE 0.1253018 0.11636026
## OTOGA 0.1619856 0.16315102
## CALJA 0.2087074 0.23122606
## MACMU 0.1904524 0.24372696
## NOMLE 0.2347050 0.32149102
## GORGO 0.2540529 0.36566381
## PANTR 0.2336738 0.34285142
## DENISOVAN 0.5903055 0.48754410
## NEANDERTHAL 0.1045043 0.23628585
## HUMAN 0.1066694 0.22554817
## MO_CMNGPDNH 0.1745729 0.21954930
## MOC_MNGPDNH 0.2358727 0.63212933
## MOCM_NGPDNH 0.2188350 0.48567944
## MOCMN_GPDNH 0.2133775 0.40029925
## MOCMNG_HNDP 0.4113882 15.96900842
## MOCMNGP_HND 0.2289946 0.35842950
## MOCMNGPD_HN 0.0923868 0.22637260
## MOCMNP_HNDG 0.1120567 0.62693088
## MOCMNDHN_PG 0.0984602 0.23519968
## MOCMNGPH_ND 0.1018453 0.21166113
## MOCMNGPN_HD 0.0161238 0.06998989
normalized_dnds_per_gene<-(apply(as.matrix(dnds_per_gene),2,function(x){
y<-x[!is.na(x)]
y[y>0]<-rank(y[y>0],ties.method="min",na.last="keep")/sum(!is.na(y[y>0]))
x[!is.na(x)]<-y
x
}))
normalized_dnds_per_gene[dnds_per_gene==0]<-0
colnames(normalized_dnds_per_gene)<-colnames(dnds_per_gene)
print(paste0("Column means/sd after normalization: "))
## [1] "Column means/sd after normalization: "
cbind(apply(as.matrix(normalized_dnds_per_gene),2,function(x){mean(x,na.rm=TRUE)}),apply(as.matrix(normalized_dnds_per_gene),2,function(x){sd(x,na.rm=TRUE)}),apply(as.matrix(normalized_dnds_per_gene),2,function(x){sum(x==min(x,na.rm=TRUE),na.rm=TRUE)}))
## [,1] [,2] [,3]
## MOUSE 0.48897649 0.2948124 198
## OTOGA 0.48026180 0.2992517 353
## CALJA 0.46373308 0.3068162 644
## MACMU 0.42912081 0.3193213 1253
## NOMLE 0.41644684 0.3228531 1461
## GORGO 0.35695197 0.3325483 2254
## PANTR 0.31400591 0.3328298 2702
## DENISOVAN 0.49456703 0.2917929 97
## NEANDERTHAL 0.13445761 0.2679131 803
## HUMAN 0.13883888 0.2709081 1606
## MO_CMNGPDNH 0.45788714 0.3092364 741
## MOC_MNGPDNH 0.38359235 0.3295723 1945
## MOCM_NGPDNH 0.32852826 0.3333506 2749
## MOCMN_GPDNH 0.33770972 0.3333556 2646
## MOCMNG_HNDP 0.20376987 0.3072317 2311
## MOCMNGP_HND 0.29878067 0.3315979 2775
## MOCMNGPD_HN 0.11509591 0.2523414 924
## MOCMNP_HNDG 0.10563380 0.2441723 505
## MOCMNDHN_PG 0.12269939 0.2590456 493
## MOCMNGPH_ND 0.13171697 0.2657447 1111
## MOCMNGPN_HD 0.03703704 0.1671900 51
#Plot and write not normalized data
boxplot(dnds_per_gene,main="Not normalized")
boxplot(dnds_per_gene,outline=FALSE,main="Not normalized (without outlier)")
png(paste0("../storage/",prefix,"dnds_distribution_plots/dnds_boxplot.png"),width=round(1200*0.45),height=300)
boxplot(dnds_per_gene,outline=FALSE,main="Not normalized (without outlier)")
dev.off();
## png
## 2
write.csv2(dnds_per_gene,paste0('../storage/',paste0("csvs/",prefix),"dnds_per_gene.csv"))
#Plot and write not normalized data with quantiles
par(mar=c(8,5,1,1))
boxplot(dnds_per_gene,main="Not normalized",las=2)
boxplot(dnds_per_gene,outline=FALSE,main="Not normalized (without outlier)",las=2)
n <- ncol(dnds_per_gene)
# width of each boxplot is 0.8
x0s <- 1:n - 0.4
x1s <- 1:n + 0.4
# these are the y-coordinates for the horizontal lines
# that you need to set to the desired values.
y0s <- apply(dnds_per_gene,2,function(x){quantile(x,0.9,na.rm=TRUE)})
# add segments
segments(x0 = x0s, x1 = x1s, y0 = y0s, col = "red")
png(paste0("../storage/",prefix,"dnds_distribution_plots/dnds_boxplot.png"),width=round(1200*0.45),height=300)
par(mar=c(8,5,1,1))
boxplot(dnds_per_gene,outline=FALSE,main="Not normalized (without outlier)",las=2)
dev.off();
## png
## 2
png(paste0("../storage/",prefix,"dnds_distribution_plots/dnds_boxplot_90percentile.png"),width=round(1200*0.45),height=300)
par(mar=c(8,5,1,1))
boxplot(dnds_per_gene,outline=FALSE,main="Not normalized (without outlier)",las=2)
segments(x0 = x0s, x1 = x1s, y0 = y0s, col = "red")
dev.off();
## png
## 2
#Plot and write normalized data
boxplot(normalized_dnds_per_gene,outline=FALSE,main="normalized (without outlier)")
boxplot(normalized_dnds_per_gene,outline=TRUE,main="normalized (with outlier)")
write.csv2(normalized_dnds_per_gene,paste0('../storage/',paste0("csvs/",prefix),"normalized_dnds_per_gene.csv"))
#Plot DNDS gene wide heatmap
normalized_dnds_per_gene_with_reordered<-normalized_dnds_per_gene[,c(1:19,21,20)]
heatmap.2(normalized_dnds_per_gene_with_reordered,
density.info="none",
trace="none",
scale="none",
dendrogram="row",
na.color="#F0F0F0",
Colv=NA,
#keysize = 0.8,
#key.par = list(mar=c(4,0.3,8,0.3),cex=6),
key.title=NA,
key.xlab="Ranked DN/DS",
cexCol=1,
srtCol= 90,
labRow = FALSE,
col=c(colorpanel(100,"#ffffb2","#feb24c","#bd0026")),
margins = c(20,5),
distfun = function(mat) {
edist <- dist(mat)
edist[which(is.na(edist))] <- max(edist, na.rm=TRUE) * 1.1
return(edist)
},
adjCol = c(1,1))
#just also save it as png or svg
#png(paste0("../storage/",prefix,"dnds_distribution_plots/dnds_heatmap.png"),width=round(3200*0.45),height=3200)
svg(paste0("../storage/",prefix,"dnds_distribution_plots/dnds_heatmap.svg"),width=9,height=20)
heatmap.2(normalized_dnds_per_gene_with_reordered,
density.info="none",
trace="none",
scale="none",
dendrogram="row",
na.color="#F0F0F0",
Colv=NA,
#keysize = 0.8,
#key.par = list(mar=c(4,0.3,8,0.3),cex=6),
key.title=NA,
key.xlab="Ranked DN/DS",
cexCol=1,
srtCol= 90,
labRow = FALSE,
col=c(colorpanel(100,"#ffffb2","#feb24c","#bd0026")),
margins = c(20,5),
distfun = function(mat) {
edist <- dist(mat)
edist[which(is.na(edist))] <- max(edist, na.rm=TRUE) * 1.1
return(edist)
},
adjCol = c(1,1))
dev.off()
## png
## 2
dnds_per_gene_with_reordered<-dnds_per_gene[,c(1:19,21,20)]
dnds_per_gene_with_reordered[dnds_per_gene_with_reordered>1.5]<-1.5
heatmap.2(dnds_per_gene_with_reordered,
density.info="none",
trace="none",
scale="none",
dendrogram="row",
na.color="#F0F0F0",
Colv=NA,
#keysize = 0.8,
#key.par = list(mar=c(4,0.3,8,0.3),cex=6),
key.title=NA,
key.xlab="Ranked DN/DS",
cexCol=1,
srtCol= 90,
labRow = FALSE,
col=c(colorpanel(100,"#ffffb2","#feb24c","#bd0026")),
margins = c(20,5),
distfun = function(mat) {
edist <- dist(mat)
edist[which(is.na(edist))] <- max(edist, na.rm=TRUE) * 1.1
return(edist)
},
adjCol = c(1,1))
#just also save it as png or svg
#png(paste0("../storage/",prefix,"dnds_distribution_plots/dnds_heatmap_raw.png"),width=round(3200*0.45),height=3200)
svg(paste0("../storage/",prefix,"dnds_distribution_plots/dnds_heatmap_raw.svg"),width=9,height=20)
heatmap.2(dnds_per_gene_with_reordered,
density.info="none",
trace="none",
scale="none",
dendrogram="row",
na.color="#F0F0F0",
Colv=NA,
#keysize = 0.8,
#key.par = list(mar=c(4,0.3,8,0.3),cex=6),
key.title=NA,
key.xlab="Ranked DN/DS",
cexCol=1,
srtCol= 90,
labRow = FALSE,
col=c(colorpanel(100,"#ffffb2","#feb24c","#bd0026")),
margins = c(20,5),
distfun = function(mat) {
edist <- dist(mat)
edist[which(is.na(edist))] <- max(edist, na.rm=TRUE) * 1.1
return(edist)
},
adjCol = c(1,1))
dev.off()
## png
## 2
#get names of genes
humanGeneNames<-unlist(mapIds(org.Hs.eg.db, paste(humanEntrezInAmba), 'SYMBOL', 'ENTREZID'))
## 'select()' returned 1:1 mapping between keys and columns
write.csv2(humanEntrezInAmba,paste0('../storage/',paste0("csvs/",prefix),"humanEntrezIDs.csv"))
#Generate genes lists for predicting functional maps WEIGHTED
addEntrezIDsHumanWeighted <- function(humanEntrezInAmba,gene_list_name_entrez_human,ratios){
setEntrezHuman<-c()
for(entrezIndex in (1:length(humanEntrezInAmba))){
humanEntrez<-humanEntrezInAmba[entrezIndex]
if(!is.na(humanEntrez)){
if(length(setEntrezHuman)<length(unique(c(setEntrezHuman,humanEntrez)))){
setEntrezHuman<-c(setEntrezHuman,humanEntrez)
geneName<-unique(humanGeneNames[humanEntrezInAmba==humanEntrez])
geneName<-geneName[!is.na(geneName)]
if(length(geneName)==0){
geneName<-paste0("entrez_",humanEntrez)
}else{
geneName<-geneName[1]
}
#humanName<-strsplit(geneName,"_")[[1]][1]
humanName<-geneName
humanName<-gsub("\\-", "_",humanName)
humanName<-gsub("\\.", "_",humanName)
humanName<-gsub("\\ ", "_",humanName)
if(!is.na(ratios[entrezIndex]) && (ratios[entrezIndex]>0)){ #exclude genes with 0 weightings (just makes the list smaller and therefore the algorithm more perfomant)
gene_list_name_entrez_human<-rbind(gene_list_name_entrez_human,c(humanName,humanEntrez,ratios[entrezIndex]))
}
}
}
}
return(gene_list_name_entrez_human)
}
gene_list_name_entrez_human<-matrix("",0,3)
set.seed(42)
for(branch in transistionsToUseForFunctionalMap){
print(paste0("branch: ",branch))
print(paste0("GeneAmount:",sum(!is.na(normalized_dnds_per_gene[,branch]) & (normalized_dnds_per_gene[,branch]>0))))
gene_list_name_entrez_human<-rbind(gene_list_name_entrez_human,c(paste0("dndsweighted_trans",branch),"",""))
gene_list_name_entrez_human<-addEntrezIDsHumanWeighted(humanEntrezInAmba,gene_list_name_entrez_human,(normalized_dnds_per_gene[,branch]))
gene_list_name_entrez_human<-rbind(gene_list_name_entrez_human,c("","",""))
}
## [1] "branch: 11"
## [1] "GeneAmount:8045"
## [1] "branch: 12"
## [1] "GeneAmount:6405"
## [1] "branch: 13"
## [1] "GeneAmount:5264"
## [1] "branch: 14"
## [1] "GeneAmount:5503"
## [1] "branch: 15"
## [1] "GeneAmount:1588"
## [1] "branch: 16"
## [1] "GeneAmount:4118"
## [1] "branch: 17"
## [1] "GeneAmount:275"
## [1] "branch: 10"
## [1] "GeneAmount:616"
dir.create("../storage/setsToRunHuman")
## Warning in dir.create("../storage/setsToRunHuman"): '..\storage\setsToRunHuman'
## already exists
write.table(gene_list_name_entrez_human,paste0(paste0("../storage/setsToRunHuman",'/',strsplit(functionalMapsHuman,"/")[[1]][length(strsplit(functionalMapsHuman,"/")[[1]])],".csv")),sep=";", row.names = FALSE, col.names = FALSE, quote = FALSE)
#Load task fmri data
taskfmri<-readMat(paste0("../storage/taskfMRI.mat"))$taskfMRI
tasknames<-unlist(readMat(paste0("../storage/taskfMRI.mat"))$tasknames)
tasknames <- gsub("tfMRI_", "", tasknames)
#Remove the contrast tasks i.e. only take the direct tasks
dif_tasks<-c(3,6,9,10,17:22,25,28,39,44:47)
tasknames<-tasknames[-dif_tasks]
taskfmri<-taskfmri[,-dif_tasks]
colnames(taskfmri)<-tasknames
#define taskgroups
colnames(taskfmri)[5:6] <- paste("LA-AVG", colnames(taskfmri)[5:6], sep = "_")
colnames(taskfmri)[c(7,9)] <- paste("MO-RF+LF", colnames(taskfmri)[c(7,9)], sep = "_")
colnames(taskfmri)[11] <- paste("MO-T", colnames(taskfmri)[11], sep = "_")
colnames(taskfmri)[c(8,10)] <- paste("MO-RH+LH", colnames(taskfmri)[c(8,10)], sep = "_")
colnames(taskfmri)[12] <- paste("MO-RH+LH", colnames(taskfmri)[12], sep = "_")
colnames(taskfmri)[16] <- paste("SO-TOM", colnames(taskfmri)[16], sep = "_")
colnames(taskfmri)[1] <- paste("EM-FACES", colnames(taskfmri)[1], sep = "_")
colnames(taskfmri)[2] <- paste("EM-SHAPE", colnames(taskfmri)[2], sep = "_")
taskfmri[,3]<-apply(taskfmri[,3:4],1,mean)
colnames(taskfmri)[4] <- paste("GA-REWARD", colnames(taskfmri)[4], sep = "_")
colnames(taskfmri)[3] <- paste("GA-AVG", "GAMBLING_AVG", sep = "_")
colnames(taskfmri)[13:14] <- paste("RE-AVG", colnames(taskfmri)[13:14], sep = "_")
colnames(taskfmri)[c(18,22,28)] <- paste("WM-FACE", colnames(taskfmri)[c(18,22,28)], sep = "_")
colnames(taskfmri)[c(17,19:21,23:27,29:30)] <- paste("WM-REST", colnames(taskfmri)[c(17,19:21,23:27,29:30)], sep = "_")
#Remove tasks that were not used in the paper
taskfmri<-taskfmri[,-15]
taskfmri<-taskfmri[,-2]
taskfmri[is.na(taskfmri)]<-0
taskfmri[abs(taskfmri)<5]<-0 #removing insiginifcant values (smaller than 5 sigmas)
heatmap.2(cor(taskfmri),cexCol=0.5,cexRow=0.5,scale="none",density.info="none", trace="none")
colorsOfBiopsyHuman<-c()
for(actRegion in 1:length(atlasRegionsHuman)){
colorsOfBiopsyHuman<-c(colorsOfBiopsyHuman,paste0("#",getAcronymByID(documentHuman$msg[[1]],atlasRegionsHuman[[actRegion]][1])$color_hex_triplet))
}
namesOfBiopsyHuman<-c()
for(actRegion in 1:length(atlasRegionsHuman)){
namesOfBiopsyHuman<-c(namesOfBiopsyHuman,paste0(getAcronymByID(documentHuman$msg[[1]],atlasRegionsHuman[[actRegion]][1])$name))
}
taskClasses<-sapply(colnames(taskfmri),function(x){sapply(strsplit(x, "_"),"[[",1)})
taskGroupfmri<-sapply(unique(taskClasses),function(x){
if(sum(x==taskClasses)>1){
rowMeans(taskfmri[,x==taskClasses])
}else{
taskfmri[,x==taskClasses]
}})
rownames(taskGroupfmri)<-namesOfBiopsyHuman
# heatmap.2(t(taskGroupfmri[order(colorsOfBiopsyHuman),]),hclustfun=function(x){hclust(x,method="ward.D2")},main="taskGroupFMRI cortical regions",scale="none",density.info="none", trace="none",ColSideColors = sort(colorsOfBiopsyHuman),Colv=NA,col=colorpanel(100, "blue","white", "red"),margins=c(6,10),cexRow=0.8,cexCol=0.1,dendrogram="row")
#
# heatmap.2(t(abs(taskGroupfmri)[order(colorsOfBiopsyHuman),]),hclustfun=function(x){hclust(x,method="ward.D2")},main="taskGroupFMRI cortical regions",scale="none",density.info="none", trace="none",ColSideColors = sort(colorsOfBiopsyHuman),Colv=NA,col=colorpanel(100, "white","orange", "red"),margins=c(6,10),cexRow=0.8,cexCol=0.1,dendrogram="row")
#Load network data
#Creating networks from literature
networks<-matrix(0,nrow=sum(atlasRegionsHuman>0),ncol=11)
colnames(networks)<-c("cortico-limbic_network","mPFC-AMY_network","mPFC-Acb_network","deafult_mode_network","salience_network","central_executive_network","sensorimotor_network","visuospatial_network","fronto-parietal","dorsal_attention","ventral_attention")
networksIDHumanList<-list()
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4035,4021,4050,4028,4268,4220,4249,9391,4212,4205,4156,4096,4111,4118,4071,4125,4278,4287,4293,4392,4140,4147,4053,4327)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4021,4050,4047,4897,4062,4327)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4021,4050,4047,4897,4062,4290)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4021,4050,4047,4897,4062,4235,4228,4242,4147,4111,4249)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4225,4272,4290,9066,9072,4540,12927,4320,4313)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4021,4028,4085,9372,4278,12920)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4021,4272,4010,4085)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4028,4085,4181,4118,4184)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4028,4096,4103)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4104,4096,4028)
networksIDHumanList[[length(networksIDHumanList)+1]]<-c(4035,4133,4103)
for(actNetwork in 1:length(networksIDHumanList)){
networks[,actNetwork]<-rep(FALSE,length(atlasRegionsHuman))
for(actSubRegionID in networksIDHumanList[[actNetwork]]){
networks[,actNetwork]<-networks[,actNetwork]|(getAtlasRegionsOfIDHuman(atlasRegionsHuman,actSubRegionID)>0)
}
}
#heatmap.2(t(abs(networks)[order(colorsOfBiopsyHuman),]),hclustfun=function(x){hclust(x,method="ward.D2")},main="literature network regions",scale="none",density.info="none", trace="none",ColSideColors = sort(colorsOfBiopsyHuman),Colv=NA,col=colorpanel(100, "white","orange", "red"),margins=c(6,10),cexRow=0.8,cexCol=0.1,dendrogram="row")
#Visualize branches to networks and tasks combined on region level
regionsIDsHumanLarge<-list()
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4035)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4021)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4050)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4028)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4268)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4220)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4249)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4212)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4205)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4156)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4096)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4111)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4118)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4071)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4125)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4140)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4147)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4053)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4047)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4897)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4062)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4235)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4228)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4225)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4085)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4010)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4181)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4184)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4278)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4287)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4293)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4392)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4327)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4290)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(9066)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(9072)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4540)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4313)
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-c(4321)
namesOfRegionsLarge<-c("IFG",
"SFG",
"MOrG",
"MFG",
"Ins",
"CgG",
"HiF",
"SOG",
"IOG",
"FuG",
"SPL",
"AnG",
"Pcu",
"PCLa",
"PCLp",
"MTG",
"ITG",
"AOrG",
"GRe",
"SRoG",
"SCG",
"CgGr",
"CgGp",
"CgGf",
"PoG",
"PrG",
"OP",
"Cun",
"Cd",
"Pu",
"GP",
"TH",
"Amg",
"Acb",
"VTA",
"SN",
"Hy",
"BSTl",
"Cl")
colorsOfRegionsLarge<-c()
combinedNetworks<-cbind(abs(taskGroupfmri),abs(networks))
combinedNetworksRegionWise<-matrix(0,0,ncol(combinedNetworks))
for(actRegionIndex in 1:length(regionsIDsHumanLarge)){
actRegionAtlas<-getAtlasRegionsOfIDHuman(atlasRegionsHuman,regionsIDsHumanLarge[actRegionIndex][1])>0
a<-getAcronymByID(documentHuman$msg[[1]],regionsIDsHumanLarge[actRegionIndex][1])
colorsOfRegionsLarge<-c(colorsOfRegionsLarge,paste0("#",a$color_hex_triplet))
combinedNetworksRegionWise<-rbind(combinedNetworksRegionWise,apply(combinedNetworks,2,function(x){
max(x[actRegionAtlas])
}))
}
rownames(combinedNetworksRegionWise)<-namesOfRegionsLarge
write.csv2(t(combinedNetworksRegionWise[order(colorsOfRegionsLarge),1:ncol(abs(taskGroupfmri))]),paste0('../storage/',paste0("csvs/"),"task_network_regions.csv"))
heatmap.2(t(combinedNetworksRegionWise[order(colorsOfRegionsLarge),1:ncol(abs(taskGroupfmri))]),hclustfun=function(x){hclust(x,method="ward.D2")},main="task network regions",scale="none",density.info="none", trace="none",ColSideColors = sort(colorsOfRegionsLarge),Colv=NA,col=colorpanel(100, "white","orange", "red"),margins=c(6,10),cexRow=0.8,cexCol=0.8,dendrogram="row")
write.csv2(t(combinedNetworksRegionWise[order(colorsOfRegionsLarge),(1:ncol(abs(taskGroupfmri)))+ncol(abs(networks))]),paste0('../storage/',paste0("csvs/"),"literuature_network_regions.csv"))
heatmap.2(t(combinedNetworksRegionWise[order(colorsOfRegionsLarge),(1:ncol(abs(taskGroupfmri)))+ncol(abs(networks))]),hclustfun=function(x){hclust(x,method="ward.D2")},main="literature network regions",scale="none",density.info="none", trace="none",ColSideColors = sort(colorsOfRegionsLarge),Colv=NA,col=colorpanel(100, "white","orange", "red"),margins=c(6,10),cexRow=0.8,cexCol=0.8,dendrogram="row")
#Create gene-wise correlation to tasks and networks for biclustering
resultsTableCorrelation<-matrix(0,length(humanEntrezInAmba),ncol(taskGroupfmri)+ncol(networks))
rownames(resultsTableCorrelation)<-humanEntrezInAmba
colnames(resultsTableCorrelation)<-c(colnames(taskGroupfmri),colnames(networks))
for(k in 1:length(humanEntrezInAmba)){
resultsTableCorrelation[k,1:ncol(taskGroupfmri)]<-apply(abs(taskGroupfmri),2,function(x){
cor(x,expressionHuman[,match(humanEntrezInAmba[k],entrezIDsOfRowsHuman)],use ="pairwise.complete.obs",method="spearman")
})
resultsTableCorrelation[k,(1:ncol(networks))+ncol(taskGroupfmri)]<-apply(abs(networks),2,function(x){
cor(x,expressionHuman[,match(humanEntrezInAmba[k],entrezIDsOfRowsHuman)],use ="pairwise.complete.obs",method="spearman")
})
}
write.csv2(resultsTableCorrelation,paste0('../storage/',paste0("csvs/",prefix),"gene_task_network_correlation_brainwide.csv"))
head(resultsTableCorrelation)
## EM-FACES GA-AVG GA-REWARD LA-AVG MO-RF+LF MO-RH+LH
## 10043 -0.12225419 -0.13932558 -0.13747772 -0.21584027 -0.09137790 -0.142723895
## 1004 -0.06171345 -0.04170594 -0.03138104 -0.07768244 0.05390198 0.004952909
## 10083 -0.12778216 -0.10381126 -0.10113888 -0.19603992 -0.06891809 -0.115020047
## 10160 0.04232188 0.07811879 0.04979313 0.20333776 -0.07031606 -0.044040218
## 1016 0.19970219 0.18264748 0.18495099 0.26711596 0.16196002 0.221450100
## 10270 0.19452623 0.24256334 0.23319998 0.39914308 0.22099737 0.283286716
## MO-T RE-AVG SO-TOM WM-REST WM-FACE
## 10043 -0.09264273 -0.112950052 -0.10688686 -0.16614041 -0.15502639
## 1004 0.01075959 -0.034766388 -0.03058994 -0.02582245 -0.03220408
## 10083 -0.08372677 -0.075708591 -0.15120910 -0.12257321 -0.10260770
## 10160 0.04798651 -0.005703194 0.08324693 0.07407935 0.04643326
## 1016 0.15469942 0.195537701 0.16291591 0.21822855 0.18827356
## 10270 0.20961120 0.188216037 0.16667897 0.29194672 0.26205380
## cortico-limbic_network mPFC-AMY_network mPFC-Acb_network
## 10043 -0.15934794 -0.06373730 -0.074460375
## 1004 0.11933947 -0.03413417 0.005385308
## 10083 -0.05207671 -0.11742923 -0.093627618
## 10160 0.28324399 0.20443041 0.160412495
## 1016 -0.07207181 0.04023532 0.071913270
## 10270 0.26749063 0.09708940 0.143436205
## deafult_mode_network salience_network central_executive_network
## 10043 -0.09128415 0.09632715 -0.079311424
## 1004 -0.16226170 0.07869912 0.175555986
## 10083 -0.07432311 -0.05460801 -0.003010283
## 10160 0.29344091 0.12010103 0.010111585
## 1016 -0.01010257 -0.08807999 0.018053726
## 10270 0.16956488 -0.19097613 0.193567471
## sensorimotor_network visuospatial_network fronto-parietal
## 10043 -0.04617075 -0.025183888 -0.05240495
## 1004 0.09820619 0.006987545 -0.01184755
## 10083 0.00477559 -0.016036088 -0.10640415
## 10160 0.06616453 0.036443746 0.08836849
## 1016 0.10784157 0.115120549 0.11087893
## 10270 0.15316757 0.113855032 0.18731775
## dorsal_attention ventral_attention
## 10043 -0.05721863 -0.057332473
## 1004 -0.00368351 0.003709607
## 10083 -0.08159040 -0.103653613
## 10160 0.08971130 0.099709028
## 1016 0.10450219 0.103227171
## 10270 0.17103961 0.186488822
gene_wise_network_correlation<-resultsTableCorrelation
gene_wise_network_correlationOrg<-gene_wise_network_correlation
gene_wise_network_correlation[,1:11]<-rank(gene_wise_network_correlation[,1:11],ties.method="max",na.last=FALSE)/length(unlist(gene_wise_network_correlation[,1:11]))
gene_wise_network_correlation[,12:22]<-rank(gene_wise_network_correlation[,12:22],ties.method="max",na.last=FALSE)/length(unlist(gene_wise_network_correlation[,12:22]))
gene_wise_network_correlation[,]<-t(apply(gene_wise_network_correlation[,],1,function(x){
rank(x,ties.method="max",na.last=FALSE)/length(x)
}))
gene_wise_network_correlation[,]<-t(apply(gene_wise_network_correlation[,],1,function(x){
x[x==min(x)]<-0
x
}))
gene_wise_network_correlation[gene_wise_network_correlationOrg<=0.1]<-0
#DNDS to network correlation
time_clusters_celltype_cor<-cor(dnds_per_gene,(gene_wise_network_correlation),method="spearman",use="pairwise.complete.obs")
rownames(time_clusters_celltype_cor)<-colnames(dnds_per_gene)
colnames(time_clusters_celltype_cor)<-colnames(gene_wise_network_correlation)
time_clusters_celltype_cor<-time_clusters_celltype_cor[,apply(time_clusters_celltype_cor,2,function(x){sum(is.na(x))})<nrow(time_clusters_celltype_cor)]
time_clusters_celltype_cor[is.na(time_clusters_celltype_cor)]<-0
heatmap.2((time_clusters_celltype_cor),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "blue","white", "red"),key.xlab = "Spearman rank correlation",dendrogram="col",Rowv=NA)
dev.copy(png,paste0("../storage/biclustering/DNDS_col_Network_col_cor.png"),width=1000,height=500)
## png
## 3
dev.off()
## png
## 2
write.csv2(time_clusters_celltype_cor,paste0("../storage/biclustering/DNDS_col_Network_col_cor.csv"))
time_clusters_celltype_cor_norm<-t(apply(time_clusters_celltype_cor,1,scale))
time_clusters_celltype_cor_norm[time_clusters_celltype_cor_norm>3]<-3
time_clusters_celltype_cor_norm[time_clusters_celltype_cor_norm<(-3)]<-(-3)
colnames(time_clusters_celltype_cor_norm)<-colnames(time_clusters_celltype_cor)
rownames(time_clusters_celltype_cor_norm)<-rownames(time_clusters_celltype_cor)
heatmap.2((time_clusters_celltype_cor_norm),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "blue","white", "red"),key.xlab = "Row z-score",dendrogram="col",Rowv=NA)
dev.copy(png,paste0("../storage/biclustering/DNDS_col_Network_col_cor_row_norm.png"),width=1000,height=500)
## png
## 3
dev.off()
## png
## 2
# Calculate functional maps
############ RUN UNTIL HERE, THEN COMPUTE p-value maps ############
############ YOU SHOULD BETTER RUN THIS ON A MACHINE ############
############ WITH >10 CORES SINCE THIS CAN TAKE A WHILE ############
#you don't need to run this, results are included
#source("predict_human_expression_synergy.R")
#Compare functional maps (p-value maps) of branches to taskfMRI
resultsTableCorrelation<-matrix(0,length(transistionsToUseForFunctionalMap),ncol(taskGroupfmri))
rownames(resultsTableCorrelation)<-transistionsToUseForFunctionalMap
colnames(resultsTableCorrelation)<-colnames(taskGroupfmri)
#load pvalue maps an correlate it with teh results table
actRow<-1
for(k in transistionsToUseForFunctionalMap){
if(file.exists(paste0(functionalMapsHuman,'/','pvalues/dndsweighted_trans',k,'_ttest_geneExpr.rds'))){
pvalsExprHuman<--log10(readRDS(paste0(functionalMapsHuman,'/','pvalues/dndsweighted_trans',k,'_ttest_geneExpr.rds')))
pvalsExprHuman[pvalsExprHuman>3]<-3
resultsTableCorrelation[actRow,]<-apply(abs(taskGroupfmri),2,function(x){
cor(x,pvalsExprHuman,use ="pairwise.complete.obs",method="spearman")
})
actRow<-actRow+1
}
}
rownames(resultsTableCorrelation)<-colnames(dnds_per_gene)[transistionsToUseForFunctionalMap]
heatmap.2(t(resultsTableCorrelation),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "blue","white", "red"),key.xlab = "Spearman rank correlation",dendrogram="row",Colv=NA)
dev.copy(png,paste0("../storage/",prefix,"branch_to_network_and_tasks_correlation_dndsweighted/additional_plots/branch_to_task_rank_correlation.png"))
## png
## 3
dev.off()
## png
## 2
heatmap.2(apply(resultsTableCorrelation,1,function(x){
rank(x,ties.method="max",na.last=FALSE)/length(x)
}),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "white","orange", "red"),key.xlab = "Column rank normalization",dendrogram="row",Colv=NA)
dev.copy(png,paste0("../storage/",prefix,"branch_to_network_and_tasks_correlation_dndsweighted/additional_plots/branch_to_task_rank_correlation_column_normalized.png"))
## png
## 3
dev.off()
## png
## 2
#Compare branches to networks
resultsTableCorrelation<-matrix(0,length(transistionsToUseForFunctionalMap),ncol(networks))
rownames(resultsTableCorrelation)<-transistionsToUseForFunctionalMap
colnames(resultsTableCorrelation)<-colnames(networks)
actRow<-1
for(k in transistionsToUseForFunctionalMap){
if(file.exists(paste0(functionalMapsHuman,'/','pvalues/dndsweighted_trans',k,'_ttest_geneExpr.rds')) && file.exists(paste0(functionalMapsHuman,'/','pvalues/dndsweighted_trans',k,'_ttest_geneExpr.rds'))){
pvalsExprHuman<--log10(readRDS(paste0(functionalMapsHuman,'/','pvalues/dndsweighted_trans',k,'_ttest_geneExpr.rds')))
pvalsExprHuman[pvalsExprHuman>3]<-3
resultsTableCorrelation[actRow,]<-apply(abs(networks),2,function(x){
cor(x,pvalsExprHuman,use ="pairwise.complete.obs",method="spearman")
})
}
actRow<-actRow+1
}
rownames(resultsTableCorrelation)<-colnames(dnds_per_gene)[transistionsToUseForFunctionalMap]
heatmap.2(t(resultsTableCorrelation),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "blue","white", "red"),key.xlab = "Spearman rank correlation",dendrogram="row",Colv=NA)
dev.copy(png,paste0("../storage/",prefix,"branch_to_network_and_tasks_correlation_dndsweighted/additional_plots/branch_to_network_rank_correlation.png"))
## png
## 3
dev.off()
## png
## 2
heatmap.2(apply(resultsTableCorrelation,1,function(x){
rank(x,ties.method="max",na.last=FALSE)/length(x)
}),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "white","orange", "red"),key.xlab = "Column rank normalization",dendrogram="row",Colv=NA)
dev.copy(png,paste0("../storage/",prefix,"branch_to_network_and_tasks_correlation_dndsweighted/additional_plots/branch_to_network_rank_correlation_column_normalized.png"))
## png
## 3
dev.off()
## png
## 2
#Compare branches to networks and tasks combined
resultsTableCorrelation<-matrix(0,length(transistionsToUseForFunctionalMap),ncol(networks)+ncol(taskGroupfmri))
rownames(resultsTableCorrelation)<-transistionsToUseForFunctionalMap
colnames(resultsTableCorrelation)<-c(colnames(networks),colnames(taskGroupfmri))
actRow<-1
for(k in transistionsToUseForFunctionalMap){
if(file.exists(paste0(functionalMapsHuman,'/','pvalues/dndsweighted_trans',k,'_ttest_geneExpr.rds'))){
pvalsExprHuman<--log10(readRDS(paste0(functionalMapsHuman,'/','pvalues/dndsweighted_trans',k,'_ttest_geneExpr.rds')))
pvalsExprHuman[pvalsExprHuman>3]<-3
resultsTableCorrelation[actRow,]<-c(apply(abs(networks),2,function(x){
cor(x,pvalsExprHuman,use ="pairwise.complete.obs",method="spearman")
}),
apply(abs(taskGroupfmri),2,function(x){
cor(x,pvalsExprHuman,use ="pairwise.complete.obs",method="spearman")
})
)
}
actRow<-actRow+1
}
resultsTableCorrelation[,1:ncol(networks)]<-rank(resultsTableCorrelation[,1:ncol(networks)],ties.method="max",na.last=FALSE)/length(unlist(resultsTableCorrelation[,1:ncol(networks)]))
resultsTableCorrelation[,(1:ncol(taskGroupfmri))+ncol(networks)]<-rank(resultsTableCorrelation[,(1:ncol(taskGroupfmri))+ncol(networks)],ties.method="max",na.last=FALSE)/length(unlist(resultsTableCorrelation[,(1:ncol(taskGroupfmri))+ncol(networks)]))
rownames(resultsTableCorrelation)<-colnames(dnds_per_gene)[transistionsToUseForFunctionalMap]
heatmap.2(t(resultsTableCorrelation),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "blue","white", "red"),key.xlab = "Block-ranks",dendrogram="row",Colv=NA)
dev.copy(png,paste0("../storage/",prefix,"branch_to_network_and_tasks_correlation_dndsweighted/branch_to_combined_rank_correlation.png"))
## png
## 3
dev.off()
## png
## 2
heatmap.2(t(apply(resultsTableCorrelation,2,function(x){
rank(x,ties.method="max",na.last=FALSE)/length(x)
})),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "white","orange", "red"),key.xlab = "Row rank normalization",dendrogram="row",Colv=NA)
dev.copy(png,paste0("../storage/",prefix,"branch_to_network_and_tasks_correlation_dndsweighted/additional_plots/branch_to_combined_rank_correlation_row_normalized.png"))
## png
## 3
dev.off()
## png
## 2
heatmap.2(apply(resultsTableCorrelation,1,function(x){
rank(x,ties.method="max",na.last=FALSE)/length(x)
}),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "white","orange", "red"),key.xlab = "Column rank normalization",dendrogram="row",Colv=NA)
dev.copy(png,paste0("../storage/",prefix,"branch_to_network_and_tasks_correlation_dndsweighted/branch_to_combined_rank_correlation_column_normalized.png"))
## png
## 3
dev.off()
## png
## 2
#Set large brain regions
regionNamesList<-c("myelencephalon","cerebellum","pons","pretectal region","hypothalamus","subthalamus","thalamus","habenular nuclei","paraventricular nuclei of thalamus","pineal gland","basal forebrain","basal ganglia","claustrum","fusiform gyrus","Heschl's gyrus","inferior temporal gyrus","middle temporal gyrus","planum polare","planum temporale","superior temporal gyrus","temporal pole","transverse gyri","inferior parietal lobule","paracentral lobule, posterior part","postcentral gyrus","precuneus","superior parietal lobule","cuneus","inferior occipital gyrus","lingual gyrus","occipital pole","occipito-temporal gyrus","superior occipital gyrus","cingulate gyrus","hippocampal formation","parahippocampal gyrus","piriform cortex","anterior orbital gyrus","frontal operculum","frontal pole","gyrus rectus","inferior frontal gyrus","inferior rostral gyrus","lateral orbital gyrus","medial orbital gyrus","middle frontal gyrus","paracentral lobule, anterior part","parolfactory gyri","posterior orbital gyrus","precentral gyrus","superior frontal gyrus","superior rostral gyrus","insula","amygdalohippocampal transition zone","basolateral nucleus","basomedial nucleus","central nucleus","cortico-medial group","lateral nucleus","central gray substance of midbrain","Edinger-Westphal nucleus","interstitial nucleus of Cajal","midbrain raphe nuclei","midbrain reticular formation","nucleus of Darkschewitsch","oculomotor nuclear complex","red nucleus","substantia nigra","trochlear nucleus","ventral tegmental area","inferior colliculus","superior colliculus")
regionsIDsHumanLarge<-list()
hasBiopsySites<-c()
for(actName in regionNamesList){
idOfName<-getIDbyName(documentHuman$msg[[1]],actName)$id
regionsIDsHumanLarge[[length(regionsIDsHumanLarge)+1]]<-idOfName
amountOfBiopsySites<-sum(getAtlasRegionsOfIDHuman(atlasRegionsHuman,idOfName)>0)
if(amountOfBiopsySites==0){
print(paste(actName,"with id",idOfName,"has no biopsySites!"))
}
hasBiopsySites<-c(hasBiopsySites,amountOfBiopsySites>0)
}
print(paste("Length of regionNamesList:",length(regionNamesList)," should be equal to length of regionsIDsHumanLarge:",length(regionsIDsHumanLarge)))
## [1] "Length of regionNamesList: 72 should be equal to length of regionsIDsHumanLarge: 72"
regionNamesList<-regionNamesList[hasBiopsySites]
regionsIDsHumanLarge<-regionsIDsHumanLarge[hasBiopsySites]
#Functions and precomputation for plotting human data
recursiveColorChange<-function(actNode,structureIDs,colorsOfStructures,restWhite=FALSE){
for(actChild in 1:xmlSize(actNode)){
attribs<-xmlAttrs(actNode[[actChild]])
if(!is.na(attribs["style"])){
isStructure<-attribs["structure_id"]==structureIDs
if(sum(isStructure,na.rm=TRUE)>0){
#print(attribs["style"])
if( !(substr(attribs["style"],nchar(attribs["style"])-6,nchar(attribs["style"]))=="#f2f1f0")){
if(!(substr(attribs["style"],nchar(attribs["style"])-6,nchar(attribs["style"]))=="#231f20")){
if(!(substr(attribs["style"],nchar(attribs["style"])-6,nchar(attribs["style"]))=="#6d6e70")){
if(!substr(attribs["style"],nchar(attribs["style"])-6,nchar(attribs["style"]))=="#872601"){
if(is.na(colorsOfStructures[isStructure])){
if(restWhite) attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#ffffff")
else attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#cdcdcd")
}else attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),colorsOfStructures[isStructure])
}else{
if(restWhite) attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#ffffff")
else attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#231f20")
}
}else{
if(restWhite) attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#ffffff")
else attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#6d6e70")
}
}else{
if(restWhite) attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#ffffff")
else attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#231f20")
}
}else{
if(restWhite) attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#ffffff")
else attribs["style"]<-paste0(substr(attribs["style"],0,nchar(attribs["style"])-7),"#ffffff")
}
xmlAttrs(actNode[[actChild]])<-attribs
}else{
print("BUG")
}
}
if(xmlSize(actNode[[actChild]])>0){
recursiveColorChange(actNode[[actChild]],structureIDs,colorsOfStructures,restWhite)
}
}
}
recursiveGetAllStructures<-function(actNode){
structureIDs<-c()
for(actChild in 1:xmlSize(actNode)){
attribs<-xmlAttrs(actNode[[actChild]])
if(!is.na(attribs["style"])){
structureIDs<-c(attribs["structure_id"],structureIDs)
}
if(xmlSize(actNode[[actChild]])>0){
structureIDs<-unique(c(structureIDs,recursiveGetAllStructures(actNode[[actChild]])))
}
}
return(structureIDs)
}
getSmallestParentAtlasRegions <- function(atlasRegionsHuman,ID){
if(is.null(ID)){
return(rep(FALSE,length(atlasRegionsHuman)))
}
newAtlasRegions<-getAtlasRegionsOfIDHuman(atlasRegionsHuman,ID)
if(sum(newAtlasRegions>0)==0){
newAtlasRegions<-getSmallestParentAtlasRegions(atlasRegionsHuman,getAcronymByID(documentHuman$msg[[1]],ID)$parent_structure_id)
}
return(newAtlasRegions)
}
uniqueRegionsHuman<-c()
for (x in 1:104){
fileName <- paste0("../storage/slice_svg/",x,".svg")
if(!file.exists(fileName)){
if(!dir.exists("../storage/slice_svg")){
dir.create("../storage/slice_svg")
}
if(!file.exists("../storage/slice_svg/slices.xml")){
download.file(paste0("http://api.brain-map.org/api/v2/data/query.xml?criteria=model::AtlasImage,rma::criteria,[annotated$eqtrue],alternate_images[image_type$eq%27Atlas+-+Human%27],rma::options[order$eq%27sub_images.section_number%27][num_rows$eqall]"),"../storage/slice_svg/slices.xml")
}
slicesXML=xmlParse("../storage/slice_svg/slices.xml")
slicesXMLTop = xmlRoot(slicesXML)
slicesID<-xmlValue(slicesXMLTop[['atlas-images']][[x]][["id"]][1]$text)
download.file(paste0("http://api.brain-map.org/api/v2/svg_download/",slicesID,"?groups=265297119,266932194,266932196,266932197"),fileName)
}
xmlfile=xmlParse(fileName)
xmltop = xmlRoot(xmlfile)
uniqueRegionsHuman<-unique(c(uniqueRegionsHuman,recursiveGetAllStructures(xmltop[['g']])))
}
uniqueRegionsHumanAtlas<-list()
for(actRegion in 1:length(uniqueRegionsHuman)){
uniqueRegionsHumanAtlas[[actRegion]]<-getSmallestParentAtlasRegions(atlasRegionsHuman,uniqueRegionsHuman[actRegion])
}
Mode <- function(x) {
ux <- unique(x[!is.na(x) & x>0])
if(length(ux)==0){
return(0)
}
ux[which.max(tabulate(match(x[!is.na(x) & x>0], ux)))]
}
plot_mri_style_image_human<-function(path_to_file,plotVals,maxVal,mip,isAlreadyInUniqueRegionsFormat=FALSE,discreteColorScale=FALSE,restWhite=FALSE){
#volume that will be plotted
m<-array(0, length(uniqueRegionsHuman))
plotVals[plotVals>maxVal]<-maxVal
if(isAlreadyInUniqueRegionsFormat==FALSE){
for(actRegion in 1:length(uniqueRegionsHuman)){
if(mip==FALSE)
m[actRegion]<-Mode(plotVals[uniqueRegionsHumanAtlas[[actRegion]]])
else m[actRegion]<-mean(plotVals[uniqueRegionsHumanAtlas[[actRegion]]],na.rm=TRUE)
}
}else{
m<-plotVals
}
#m[is.nan(m)]<-0
#m[is.na(m)]<-0
if(mip){
m[m<0]<-0
}else{
m[m<1]<-NA
}
if(!is.null(path_to_file)){
png(paste0(path_to_file),width=2000,height=2500)
}
par(mfrow=c(4,3),mar = rep(0.1, 4))
rb<-c()
if(discreteColorScale==TRUE){
rb<-getColorScale(maxVal)
}else{
m<-round(m*100)
rb<-getColorScale(maxVal*100+1)
}
for (x in 1:11){
fileName <- paste0("../storage/slice_svg/",x*9,".svg")
xmlfile=xmlParse(fileName)
xmltop = xmlRoot(xmlfile)
if(discreteColorScale==TRUE){
recursiveColorChange(xmltop[['g']],uniqueRegionsHuman,rb[m],restWhite=restWhite)
}else{
recursiveColorChange(xmltop[['g']],uniqueRegionsHuman,rb[m+1],restWhite=restWhite)
}
brainSVG <- rsvg(charToRaw(as(xmlfile, "character")),width=1500)
plot(c(0,0), type="n", axes=F, xlab="", ylab="")
rasterImage(brainSVG,1,-1,2,1)
}
image.scale <- function(z, zlim, col = heat.colors(12),
breaks, horiz=TRUE, ylim=NULL, xlim=NULL, ...){
if(!missing(breaks)){
if(length(breaks) != (length(col)+1)){stop("must have one more break than colour")}
}
if(missing(breaks) & !missing(zlim)){
breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1))
}
if(missing(breaks) & missing(zlim)){
zlim <- range(z, na.rm=TRUE)
zlim[2] <- zlim[2]+c(zlim[2]-zlim[1])*(1E-3)#adds a bit to the range in both directions
zlim[1] <- zlim[1]-c(zlim[2]-zlim[1])*(1E-3)
breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1))
}
poly <- vector(mode="list", length(col))
for(i in seq(poly)){
poly[[i]] <- c(breaks[i], breaks[i+1], breaks[i+1], breaks[i])
}
xaxt <- ifelse(horiz, "s", "n")
yaxt <- ifelse(horiz, "n", "s")
if(horiz){YLIM<-c(0,1); XLIM<-range(breaks)}
if(!horiz){YLIM<-range(breaks); XLIM<-c(0,1)}
if(missing(xlim)) xlim=XLIM
if(missing(ylim)) ylim=YLIM
plot(1,1,t="n",ylim=ylim, xlim=xlim, xaxt=xaxt, yaxt=yaxt, xaxs="i", yaxs="i",cex.axis=6, ...)
for(i in seq(poly)){
if(horiz){
polygon(poly[[i]], c(0,0,1,1), col=col[i], border=NA)
}
if(!horiz){
polygon(c(0,0,1,1), poly[[i]], col=col[i], border=NA)
}
}
}
par(mgp=c(5,6,0))
if(!is.null(path_to_file)){
par(mar=c(15,5,15,5))
}else{
par(mar=c(3,1,3,1))
}
if(mip){
if(discreteColorScale==TRUE){
image.scale(c(0,((maxVal))),col=getColorScale((maxVal)))
}else{
image.scale(c(0,((maxVal))),col=getColorScale((maxVal*100)))
}
}else{
if(discreteColorScale==TRUE){
image.scale(c(1,((maxVal))),col=getColorScale((maxVal)))
}else{
image.scale(c(1,((maxVal))),col=getColorScale((maxVal*100)))
}
}
if(!is.null(path_to_file)){
dev.off()
}
}
plot_mri_style_image_human_all_slices<-function(path_to_file,plotVals,maxVal,mip,isAlreadyInUniqueRegionsFormat=FALSE,discreteColorScale=FALSE,restWhite=FALSE){
#volume that will be plotted
m<-array(0, length(uniqueRegionsHuman))
plotVals[plotVals>maxVal]<-maxVal
if(isAlreadyInUniqueRegionsFormat==FALSE){
for(actRegion in 1:length(uniqueRegionsHuman)){
if(mip==FALSE)
m[actRegion]<-Mode(plotVals[uniqueRegionsHumanAtlas[[actRegion]]])
else m[actRegion]<-mean(plotVals[uniqueRegionsHumanAtlas[[actRegion]]],na.rm=TRUE)
}
}else{
m<-plotVals
}
#m[is.nan(m)]<-0
#m[is.na(m)]<-0
if(mip){
m[m<0]<-0
}else{
m[m<1]<-NA
}
png(paste0(path_to_file),width=10000,height=4000, type="cairo")
par(mfrow=c(6,17),mar = rep(0.1, 4))
rb<-c()
if(discreteColorScale==TRUE){
rb<-getColorScale(maxVal)
}else{
m<-round(m*100)
rb<-getColorScale(maxVal*100+1)
}
for (x in 1:101){
fileName <- paste0("../storage/slice_svg/",x,".svg")
xmlfile=xmlParse(fileName)
xmltop = xmlRoot(xmlfile)
if(discreteColorScale==TRUE){
recursiveColorChange(xmltop[['g']],uniqueRegionsHuman,rb[m],restWhite=restWhite)
}else{
recursiveColorChange(xmltop[['g']],uniqueRegionsHuman,rb[m+1],restWhite=restWhite)
}
brainSVG <- rsvg(charToRaw(as(xmlfile, "character")),width=1000)
plot(c(0,0), type="n", axes=F, xlab="", ylab="")
rasterImage(brainSVG,1,-1,2,1)
}
image.scale <- function(z, zlim, col = heat.colors(12),
breaks, horiz=TRUE, ylim=NULL, xlim=NULL, ...){
if(!missing(breaks)){
if(length(breaks) != (length(col)+1)){stop("must have one more break than colour")}
}
if(missing(breaks) & !missing(zlim)){
breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1))
}
if(missing(breaks) & missing(zlim)){
zlim <- range(z, na.rm=TRUE)
zlim[2] <- zlim[2]+c(zlim[2]-zlim[1])*(1E-3)#adds a bit to the range in both directions
zlim[1] <- zlim[1]-c(zlim[2]-zlim[1])*(1E-3)
breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1))
}
poly <- vector(mode="list", length(col))
for(i in seq(poly)){
poly[[i]] <- c(breaks[i], breaks[i+1], breaks[i+1], breaks[i])
}
xaxt <- ifelse(horiz, "s", "n")
yaxt <- ifelse(horiz, "n", "s")
if(horiz){YLIM<-c(0,1); XLIM<-range(breaks)}
if(!horiz){YLIM<-range(breaks); XLIM<-c(0,1)}
if(missing(xlim)) xlim=XLIM
if(missing(ylim)) ylim=YLIM
plot(1,1,t="n",ylim=ylim, xlim=xlim, xaxt=xaxt, yaxt=yaxt, xaxs="i", yaxs="i",cex.axis=6, ...)
for(i in seq(poly)){
if(horiz){
polygon(poly[[i]], c(0,0,1,1), col=col[i], border=NA)
}
if(!horiz){
polygon(c(0,0,1,1), poly[[i]], col=col[i], border=NA)
}
}
}
par(mgp=c(3,6,0))
par(mar=c(15,5,15,5))
if(mip){
if(discreteColorScale==TRUE){
image.scale(c(0,((maxVal))),col=getColorScale((maxVal)))
}else{
image.scale(c(0,((maxVal))),col=getColorScale((maxVal*100)))
}
}else{
if(discreteColorScale==TRUE){
image.scale(c(1,((maxVal))),col=getColorScale((maxVal)))
}else{
image.scale(c(1,((maxVal))),col=getColorScale((maxVal*100)))
}
}
dev.off()
}
plot_mri_style_image_human_all_slices_and_significance<-function(path_to_file,plotPvals,fdr005,fdr01,justDrawOriginalColors=FALSE){
getColorScale<-function(amountOfColors,threshold){
greyColor<-colorpanel(amountOfColors,"#f9f9f9","#a8a8a8","#7a7a7a")
heatColor<-colorpanel(amountOfColors,"#ffffb2","#feb24c","#bd0026")
if(is.na(threshold))
return(greyColor)
if(threshold<=1 || threshold==Inf)
return(greyColor)
return(c(greyColor[1:floor(threshold)],heatColor[(floor(threshold)+1):amountOfColors]))
}
#volume that will be plotted
m<-array(0, length(uniqueRegionsHuman))
plotPvals[is.na(plotPvals)]<-1
for(actRegion in 1:length(uniqueRegionsHuman)){
m[actRegion]<-quantile(-log10(plotPvals[uniqueRegionsHumanAtlas[[actRegion]]]),0.9,na.rm=TRUE)
}
m[is.nan(m)]<-0
m[is.na(m)]<-0
m[m<0]<-0
minimumPval<--log10(10^(-8))
m[m>minimumPval]<-minimumPval
png(paste0(path_to_file),width=10000,height=4000, type="cairo")
par(mfrow=c(6,17),mar = rep(0.1, 4))
m<-round(m*1000)
rb<-getColorScale(minimumPval*1000+1,-log10(fdr01)*1000+1)
for (x in 1:101){
fileName <- paste0("../storage/slice_svg/",x,".svg")
xmlfile=xmlParse(fileName)
xmltop = xmlRoot(xmlfile)
if(justDrawOriginalColors==FALSE){
recursiveColorChange(xmltop[['g']],uniqueRegionsHuman,rb[m+1])
}
brainSVG <- rsvg(charToRaw(as(xmlfile, "character")),width=1000)
plot(c(0,0), type="n", axes=F, xlab="", ylab="")
rasterImage(brainSVG,1,-1,2,1)
}
image.scale <- function(z, zlim, col = heat.colors(12),
breaks, horiz=TRUE, ylim=NULL, xlim=NULL, ...){
if(!missing(breaks)){
if(length(breaks) != (length(col)+1)){stop("must have one more break than colour")}
}
if(missing(breaks) & !missing(zlim)){
breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1))
}
if(missing(breaks) & missing(zlim)){
zlim <- range(z, na.rm=TRUE)
zlim[2] <- zlim[2]+c(zlim[2]-zlim[1])*(1E-3)#adds a bit to the range in both directions
zlim[1] <- zlim[1]-c(zlim[2]-zlim[1])*(1E-3)
breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1))
}
poly <- vector(mode="list", length(col))
for(i in seq(poly)){
poly[[i]] <- c(breaks[i], breaks[i+1], breaks[i+1], breaks[i])
}
xaxt <- ifelse(horiz, "s", "n")
yaxt <- ifelse(horiz, "n", "s")
if(horiz){YLIM<-c(0,1); XLIM<-range(breaks)}
if(!horiz){YLIM<-range(breaks); XLIM<-c(0,1)}
if(missing(xlim)) xlim=XLIM
if(missing(ylim)) ylim=YLIM
plot(1,1,t="n",ylim=ylim, xlim=xlim, xaxt=xaxt, yaxt=yaxt, xaxs="i", yaxs="i",cex.axis=8, ...)
for(i in seq(poly)){
if(horiz){
polygon(poly[[i]], c(0,0,1,1), col=col[i], border=NA)
}
if(!horiz){
polygon(c(0,0,1,1), poly[[i]], col=col[i], border=NA)
}
}
if(!(fdr005==0)){
lines(c(-log10(fdr005),-log10(fdr005)),c(0,1),lty=3,lwd=5)
}
if(!(fdr01==0)){
lines(c(-log10(fdr01),-log10(fdr01)),c(0,1),lwd=5)
}
}
par(mgp=c(3,6,0))
par(mar=c(15,5,15,5))
image.scale(c(0,minimumPval),col=rb)
dev.off()
}
#Calcualte evolutionary age of every biopsy site in human brain
pvalueMapsRessource<-matrix(0, length(transistionsToUseForFunctionalMap),sum(atlasRegionsHuman>0))
humanFunctionalMapsLogPvalueTop<-matrix(0, length(transistionsToUseForFunctionalMap),sum(atlasRegionsHuman>0))
humanFunctionalMapsUniqueMapIsSignificant<-matrix(FALSE, length(transistionsToUseForFunctionalMap),length(uniqueRegionsHumanAtlas))
plot_mri_style_image_human_all_slices_and_significance(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/evoHumanHistory_region_colors.png"),rep(0,ncol(humanFunctionalMapsLogPvalueTop)),0,0,justDrawOriginalColors=TRUE)
## png
## 2
branch<-1
for(k in transistionsToUseForFunctionalMap){
print(paste0("branch: ",branch))
humanFunctionalMapsLogPvalueTop[branch,]<-readRDS(paste0(functionalMapsHuman,'/','pvalues/',"dndsweighted_trans",k,'_ttest_geneExpr.rds'))
sigThresh<-(-log10(max(humanFunctionalMapsLogPvalueTop[branch,][p.adjust(humanFunctionalMapsLogPvalueTop[branch,],method="BH")<=0.1])))
for(actRegion in 1:length(uniqueRegionsHumanAtlas)){
if(!is.na(sigThresh)){
if(quantile(-log10(humanFunctionalMapsLogPvalueTop[branch,][uniqueRegionsHumanAtlas[[actRegion]]]),0.9,na.rm=TRUE)>=sigThresh){
humanFunctionalMapsUniqueMapIsSignificant[branch,actRegion]<-TRUE
}
}
}
fdr005<-max(humanFunctionalMapsLogPvalueTop[branch,][p.adjust(humanFunctionalMapsLogPvalueTop[branch,],method="BH")<=0.05])
fdr01<-max(humanFunctionalMapsLogPvalueTop[branch,][p.adjust(humanFunctionalMapsLogPvalueTop[branch,],method="BH")<=0.1])
plot_mri_style_image_human_all_slices_and_significance(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/evoHumanHistory_",colsToUseNames[k],"_90percentile_significant_sites",transistionsToUseForFunctionalMapPostfix,".png"),humanFunctionalMapsLogPvalueTop[branch,],fdr005,fdr01)
pvalueMapsRessource[branch,]<-humanFunctionalMapsLogPvalueTop[branch,]
humanFunctionalMapsLogPvalueTop[branch,]<-(-log10(humanFunctionalMapsLogPvalueTop[branch,]))
branch<-branch+1
}
## [1] "branch: 1"
## [1] "branch: 2"
## [1] "branch: 3"
## [1] "branch: 4"
## [1] "branch: 5"
## [1] "branch: 6"
## [1] "branch: 7"
## [1] "branch: 8"
evolutionaryAgeUniqueRegionsHumanAtlasWiseTop<-unlist(sapply(uniqueRegionsHumanAtlas,function(x){
if(sum(x)>1){
which.max(apply(humanFunctionalMapsLogPvalueTop[,x],1,function(y){
mean(y,na.rm=TRUE)
}))
}else{
which.max(humanFunctionalMapsLogPvalueTop[,x])
}
}))
humanFunctionalMapsLogPvalueTop[humanFunctionalMapsLogPvalueTop>15]<-15
pvalueMapsRessource<-t(pvalueMapsRessource)
colnames(pvalueMapsRessource)<-colsToUseNames[transistionsToUseForFunctionalMap]
appendSheet<-FALSE
options(scipen = 999)
options(warn=-1)
if(file.exists( paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/pvalue_ressource",transistionsToUseForFunctionalMapPostfix,".xlsx"))){
file.remove( paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/pvalue_ressource",transistionsToUseForFunctionalMapPostfix,".xlsx"))
}
## [1] TRUE
for(i in 1:ncol(pvalueMapsRessource)){
regionPvals<-cbind(regionInfo,pvalueMapsRessource[,i])
colnames(regionPvals)[ncol(regionPvals)]<-"pvalue"
write.xlsx(regionPvals, file = paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/pvalue_ressource",transistionsToUseForFunctionalMapPostfix,".xlsx"),
sheetName = colnames(pvalueMapsRessource)[i], append = appendSheet,row.names = FALSE)
appendSheet<-TRUE
}
options(warn=0)
evolutionaryAgeBiopsyWiseTop<-unlist(apply((humanFunctionalMapsLogPvalueTop),2,function(x){return(which.max(x))}))
evolutionaryAgeBiopsyWiseTop<-sapply(uniqueRegionsHumanAtlas,function(x){Mode(evolutionaryAgeBiopsyWiseTop[x])})
evolutionaryAgeBiopsyWiseTopSign<-sapply(1:length(evolutionaryAgeBiopsyWiseTop),function(x){
if(is.na(evolutionaryAgeBiopsyWiseTop[x])){
return(NA)
}
if(humanFunctionalMapsUniqueMapIsSignificant[evolutionaryAgeBiopsyWiseTop[x],x]==FALSE){
return(NA)
}else{
return(evolutionaryAgeBiopsyWiseTop[x])
}
})
evolutionaryAgeBiopsyWiseTopHom<-unlist(apply((humanFunctionalMapsLogPvalueTop[(length(transistionsToUseForFunctionalMap)-2):length(transistionsToUseForFunctionalMap),]),2,function(x){return(which.max(x)+(length(transistionsToUseForFunctionalMap)-3))}))
evolutionaryAgeBiopsyWiseTopHom<-sapply(uniqueRegionsHumanAtlas,function(x){Mode(evolutionaryAgeBiopsyWiseTopHom[x])})
evolutionaryAgeBiopsyWiseTopHomSign<-sapply(1:length(evolutionaryAgeBiopsyWiseTopHom),function(x){
if(is.na(evolutionaryAgeBiopsyWiseTopHom[x])){
return(NA)
}
if(humanFunctionalMapsUniqueMapIsSignificant[evolutionaryAgeBiopsyWiseTopHom[x],x]==FALSE){
return(NA)
}else{
return(evolutionaryAgeBiopsyWiseTopHom[x])
}
})
evolutionaryAgeBiopsyWiseWATop<-apply((humanFunctionalMapsLogPvalueTop),2,function(x){
trans<-1:length(x)
if(sum(is.na(x))==length(x)){
return(NA)
}
y<-x
y[y<0]<-0
return(sum((trans[!is.na(y)])*y[!is.na(y)]/(sum(y,na.rm=TRUE))))
})
evolutionaryAgeBiopsyWiseWATop<-sapply(uniqueRegionsHumanAtlas,function(x){mean(evolutionaryAgeBiopsyWiseWATop[x],na.rm=TRUE)})
#Cluster Evo History By Gene Expression
dir.create(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted"))
## Warning in dir.create(paste0("../storage/",
## prefix, "evoHistoryHuman_dndsweighted")): '..
## \storage\NGT_SPLIT_TREE_9000_NAto0_evoHistoryHuman_dndsweighted' already exists
expressionHumanTopLargeRegions<-sapply(1:length(regionsIDsHumanLarge),function(actRegion){
actRegionAtlas<-getAtlasRegionsOfIDHuman(atlasRegionsHuman,regionsIDsHumanLarge[[actRegion]][1])>0
if(sum(actRegionAtlas)==1){
return(expressionHuman[actRegionAtlas,is.element(entrezIDsOfRowsHuman,humanEntrezInAmba)])
}else{
return(apply(expressionHuman[actRegionAtlas,is.element(entrezIDsOfRowsHuman,humanEntrezInAmba)],2,function(x){
mean(x,na.rm=TRUE)
}))
}
})
colnames(expressionHumanTopLargeRegions)<-1:ncol(expressionHumanTopLargeRegions)
for(actRegion in 1:length(regionsIDsHumanLarge)){
colnames(expressionHumanTopLargeRegions)[actRegion]<-getAcronymByID(documentHuman$msg[[1]],regionsIDsHumanLarge[[actRegion]])$acronym
}
cordist<-as.dist((1-cor(expressionHumanTopLargeRegions))/2)
cordist[is.na(cordist)]<-0
hc <- hclust(cordist,method="ward.D2")
#hc <- hclust(dist(t(expressionHumanTopLargeRegions)),method="ward.D2")
clustersExpression<-cutree(hc, k = 8)
hist(clustersExpression)
plot(hc)
dev.copy(png,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/gene_expression_clustering_dendrogram",transistionsToUseForFunctionalMapPostfix,".png"),width=900,height=400)
## png
## 3
dev.off()
## png
## 2
clustersExpressionCSV<-paste("Cluster",clustersExpression)
names(clustersExpressionCSV)<-colnames(expressionHumanTopLargeRegions)
sortedclustersExpressionCSV<-sort(clustersExpressionCSV)
write.csv2(sortedclustersExpressionCSV,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/expression_regions_cluster",transistionsToUseForFunctionalMapPostfix,".csv"))
clustersExpressionBiopsyLevel<-rep(NA,length(atlasRegionsHuman))
for(actRegion in 1:length(regionsIDsHumanLarge)){
actRegionAtlas<-getAtlasRegionsOfIDHuman(atlasRegionsHuman,regionsIDsHumanLarge[[actRegion]][1])
clustersExpressionBiopsyLevel[actRegionAtlas]<-clustersExpression[actRegion]
}
clusterProfiles<-matrix(0,dim(expressionHumanTopLargeRegions)[1],max(clustersExpression,na.rm=TRUE))
for(timepoint in 1:dim(expressionHumanTopLargeRegions)[1]){
for(cluster in unique(clustersExpression)){
clusterProfiles[timepoint,cluster]<-mean(expressionHumanTopLargeRegions[timepoint,clustersExpression==cluster],rm.na=TRUE)
}
}
library(RColorBrewer)
getColorScale<-function(amountOfColors){
#return(gsub('.{2}$', '',rainbow(amountOfColors)))
return(brewer.pal(n = 8, name = "Set2"))
}
rb<-getColorScale(max(clustersExpression,na.rm=TRUE))
plot_mri_style_image_human(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/expression_clustering",transistionsToUseForFunctionalMapPostfix,".png"),clustersExpressionBiopsyLevel,max(clustersExpressionBiopsyLevel,na.rm=TRUE),mip=FALSE,isAlreadyInUniqueRegionsFormat=FALSE,discreteColorScale=TRUE,restWhite=FALSE)
## png
## 2
plot_mri_style_image_human_all_slices(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/expression_clustering",transistionsToUseForFunctionalMapPostfix,".png"),clustersExpressionBiopsyLevel,max(clustersExpressionBiopsyLevel,na.rm=TRUE),mip=FALSE,isAlreadyInUniqueRegionsFormat=FALSE,discreteColorScale=TRUE,restWhite=FALSE)
## png
## 2
#Cluster Evo History By Functional Maps
humanFunctionalMapsLogPvalueTopLargeRegions<-sapply(1:length(regionsIDsHumanLarge),function(actRegion){
actRegionAtlas<-getAtlasRegionsOfIDHuman(atlasRegionsHuman,regionsIDsHumanLarge[[actRegion]][1])>0
if(sum(actRegionAtlas)==1){
return(humanFunctionalMapsLogPvalueTop[,actRegionAtlas])
}else{
return(apply(humanFunctionalMapsLogPvalueTop[,actRegionAtlas],1,function(x){
mean(x,na.rm=TRUE)
}))
}
})
humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA<-humanFunctionalMapsLogPvalueTopLargeRegions
humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA[is.na(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)]<-0
humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA<-t(apply(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA,1,rank)-1)
humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA<-humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA/apply(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA,1,max)
humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA[humanFunctionalMapsLogPvalueTopLargeRegions==0 | is.na(humanFunctionalMapsLogPvalueTopLargeRegions)]<-0
colnames(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)<-1:ncol(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)
for(actRegion in 1:length(regionsIDsHumanLarge)){
colnames(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)[actRegion]<-getAcronymByID(documentHuman$msg[[1]],regionsIDsHumanLarge[[actRegion]])$acronym
}
write.csv2(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_region_pvals",transistionsToUseForFunctionalMapPostfix,".csv"))
cordist<-as.dist((1-cor(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA))/2)
cordist[is.na(cordist)]<-0
hc <- hclust(dist(t(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)),method="ward.D2")
clustersEvoHist<-cutree(hc, k = 8)
hist(clustersEvoHist)
plot(hc)
dev.copy(png,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_dendrogram",transistionsToUseForFunctionalMapPostfix,".png"),width=900,height=400)
## png
## 3
dev.off()
## png
## 2
clustersEvoHistCSV<-paste("Cluster",clustersEvoHist)
names(clustersEvoHistCSV)<-colnames(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)
sortedclustersExpressionCSV<-sort(clustersEvoHistCSV)
write.csv2(sortedclustersExpressionCSV,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_cluster",transistionsToUseForFunctionalMapPostfix,".csv"))
clustersEvoHistBipsyLevel<-rep(NA,length(atlasRegionsHuman))
for(actRegion in 1:length(regionsIDsHumanLarge)){
actRegionAtlas<-getAtlasRegionsOfIDHuman(atlasRegionsHuman,regionsIDsHumanLarge[[actRegion]][1])
clustersEvoHistBipsyLevel[actRegionAtlas]<-clustersEvoHist[actRegion]
}
clusterProfiles<-matrix(0,dim(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)[1],max(clustersEvoHist,na.rm=TRUE))
for(timepoint in 1:dim(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)[1]){
for(cluster in unique(clustersEvoHist)){
clusterProfiles[timepoint,cluster]<-mean(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA[timepoint,clustersEvoHist==cluster],rm.na=TRUE)
}
}
library(RColorBrewer)
getColorScale<-function(amountOfColors){
#return(gsub('.{2}$', '',rainbow(amountOfColors)))
return(brewer.pal(n = 8, name = "Set2"))
}
rb<-getColorScale(max(clustersEvoHist,na.rm=TRUE))
matplot(1:nrow(clusterProfiles), clusterProfiles, lty=1,type='l', xlab='DNDS column', ylab='rank normalized mean p-value (selection pressure)', col=rb[1:max(clustersEvoHist,na.rm=TRUE)],ylim=c(0,1.1),lwd=2,xaxt="n")
axis(side=1,at=1:nrow(clusterProfiles),labels=transistionsToUseForFunctionalMap)
rownames(clusterProfiles)<-transistionsToUseForFunctionalMap
colnames(clusterProfiles)<-paste("Cl",1:ncol(clusterProfiles))
write.csv2(clusterProfiles,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_timeline",transistionsToUseForFunctionalMapPostfix,".csv"))
legend('topright', inset=.05, legend=paste("Cl",1:ncol(clusterProfiles)),
lty=1, horiz=TRUE, col=rb[1:max(clustersEvoHist,na.rm=TRUE)])
dev.copy(png,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_timeline",transistionsToUseForFunctionalMapPostfix,".png"),width=900,height=400)
## png
## 3
dev.off()
## png
## 2
plot_mri_style_image_human(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering",transistionsToUseForFunctionalMapPostfix,".png"),clustersEvoHistBipsyLevel,max(clustersEvoHistBipsyLevel,na.rm=TRUE),mip=FALSE,isAlreadyInUniqueRegionsFormat=FALSE,discreteColorScale=TRUE,restWhite=FALSE)
## png
## 2
plot_mri_style_image_human_all_slices(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/evoHumanHistory_clustering",transistionsToUseForFunctionalMapPostfix,".png"),clustersEvoHistBipsyLevel,max(clustersEvoHistBipsyLevel,na.rm=TRUE),mip=FALSE,isAlreadyInUniqueRegionsFormat=FALSE,discreteColorScale=TRUE,restWhite=FALSE)
## png
## 2
regionCorrelation<-cor(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)
heatmap.2(regionCorrelation,scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "blue","white", "red"),key.xlab = "Correlation",Rowv=order.dendrogram(as.dendrogram(hc)),Colv=order.dendrogram(as.dendrogram(hc)),hclustfun=function(x){hc})
dev.copy(png,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_correlation",transistionsToUseForFunctionalMapPostfix,".png"))
## png
## 3
dev.off()
## png
## 2
heatmap.2(cor(expressionHumanTopLargeRegions),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "blue","white", "red"),key.xlab = "Correlation",Rowv=order.dendrogram(as.dendrogram(hc)),Colv=order.dendrogram(as.dendrogram(hc)),hclustfun=function(x){hc})
dev.copy(png,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_expression_correlation",transistionsToUseForFunctionalMapPostfix,".png"))
## png
## 3
dev.off()
## png
## 2
heatmap.2(regionCorrelation-cor(expressionHumanTopLargeRegions),scale="none",density.info="none",cexRow=0.8,cexCol=0.8, trace="none",col=colorpanel(100, "blue","white", "red"),key.xlab = "Correlation Distance",Rowv=order.dendrogram(as.dendrogram(hc)),Colv=order.dendrogram(as.dendrogram(hc)),hclustfun=function(x){hc})
dev.copy(png,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_expression_correlation_dist",transistionsToUseForFunctionalMapPostfix,".png"))
## png
## 3
dev.off()
## png
## 2
#Network Figure
colorsOfRegions<-sapply(1:length(regionsIDsHumanLarge),function(actRegion){
paste0("#",getAcronymByID(documentHuman$msg[[1]],regionsIDsHumanLarge[[actRegion]][1])$color_hex_triplet)
})
regionCorrelation<-cor(humanFunctionalMapsLogPvalueTopLargeRegionsWithoutNA)
write.csv2(regionCorrelation,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_correlation.csv"))
links <- matrix(0,nrow=0,ncol=3)
for(i in 1:nrow(regionCorrelation)){
for(j in i:ncol(regionCorrelation)){
if(regionCorrelation[i,j]>0 && i!=j){
links <- rbind(links,c(rownames(regionCorrelation)[i],rownames(regionCorrelation)[j],regionCorrelation[i,j]))
}
}
}
links<-data.frame(
source=links[,1],
target=links[,2],
importance=as.numeric(links[,3])*100,
stringsAsFactors=FALSE
)
nodes <- data.frame(
name=rownames(regionCorrelation)
)
# Turn it into igraph object
network <- graph_from_data_frame(d=links, vertices=nodes, directed=F)
E(network)$color<-gray.colors(100,start = 1, end = 0.1)[round(links$importance)+1]
E(network)$weights<-links$importance
V(network)$frame.color <- "white"
V(network)$label.color <- "black"
V(network)$label.font=2
V(network)$label.family="sans"
layoutingAlg<-layout_with_kk(network)
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_with_fr(network))
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_with_gem(network))
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_with_graphopt(network))
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_with_kk(network))
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_with_lgl(network))
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_with_mds(network))
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_with_dh(network))
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_with_sugiyama(network))
# plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layout_nicely(network))
plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layoutingAlg)
plot(network, vertex.color=rb[clustersEvoHist],vertex.size=8, layout=layoutingAlg)
nodecolors<-cbind(rownames(regionCorrelation),colorsOfRegions,rb[clustersEvoHist])
colnames(nodecolors)<-c("regionname","AHBA color","clustercolor")
write.csv2(nodecolors,"C:\\Users\\ganglberger\\Desktop\\nodecolors.csv")
# Make the plot
svg(filename=paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_correlation_graph",transistionsToUseForFunctionalMapPostfix,".svg"),width=20, height=20)
plot(network, vertex.color=colorsOfRegions,vertex.size=8, layout=layoutingAlg)
dev.off()
## png
## 2
svg(filename=paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_correlation_graph_cluster_cols",transistionsToUseForFunctionalMapPostfix,".svg"),width=20, height=20)
plot(network, vertex.color=rb[clustersEvoHist],vertex.size=8, layout=layoutingAlg)
dev.off()
## png
## 2
networksCombined<-cbind(taskGroupfmri,networks)
colnames(networksCombined)<-c(colnames(taskGroupfmri),colnames(networks))
networksCombined_on_network_res<-sapply(1:length(regionsIDsHumanLarge),function(actRegion){
actRegionAtlas<-getAtlasRegionsOfIDHuman(atlasRegionsHuman,regionsIDsHumanLarge[[actRegion]][1])>0
if(sum(actRegionAtlas)==1){
return(networksCombined[actRegionAtlas,])
}else{
return(apply(networksCombined[actRegionAtlas,],2,function(x){
max(x,na.rm=TRUE)
}))
}
})
networksCombined_on_network_res[networksCombined_on_network_res<0]<-0
networkscombined_on_network_res_rank<-t(apply(networksCombined_on_network_res,1,function(x){rank(x,ties.method="max",na.last="keep")/length(x)}))
networkscombined_on_network_res_rank[networksCombined_on_network_res==0]<-0
for(act_network_index in 1:ncol(networksCombined)){
svg(filename=paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_clustering_correlation_graph_",colnames(networksCombined)[act_network_index],"",transistionsToUseForFunctionalMapPostfix,".svg"),width=20, height=20)
newColors<-rb[clustersEvoHist]
newColors[networkscombined_on_network_res_rank[act_network_index,]==0]<-"#FFFFFF"
plot(network, vertex.color=newColors,vertex.size=8, layout=layoutingAlg)
dev.off()
}
#Plot evolutionary age of every biopsy site in human brain
number_colnames<-1:length(transistionsToUseForFunctionalMap)
names(number_colnames)<-colnames(dnds_per_gene)[transistionsToUseForFunctionalMap]
write.csv2(number_colnames,paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_numbers_to_columnnames",transistionsToUseForFunctionalMapPostfix,".csv"))
getColorScale<-function(amountOfColors){
return(jet.colors(amountOfColors))
}
plot_mri_style_image_human(NULL,evolutionaryAgeUniqueRegionsHumanAtlasWiseTop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=TRUE,isAlreadyInUniqueRegionsFormat=TRUE)
plot_mri_style_image_human(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_regionwise_peak",transistionsToUseForFunctionalMapPostfix,".png"),evolutionaryAgeUniqueRegionsHumanAtlasWiseTop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=TRUE,isAlreadyInUniqueRegionsFormat=TRUE)
## png
## 2
plot_mri_style_image_human(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_peak_modeRegion",transistionsToUseForFunctionalMapPostfix,".png"),evolutionaryAgeBiopsyWiseTop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=TRUE,isAlreadyInUniqueRegionsFormat=TRUE)
## png
## 2
plot_mri_style_image_human(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_peak_modeRegion",transistionsToUseForFunctionalMapPostfix,"_sign.png"),evolutionaryAgeBiopsyWiseTopSign,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=TRUE,isAlreadyInUniqueRegionsFormat=TRUE)
## png
## 2
plot_mri_style_image_human(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistoryWA_meanRegion",transistionsToUseForFunctionalMapPostfix,".png"),evolutionaryAgeBiopsyWiseWATop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=FALSE,isAlreadyInUniqueRegionsFormat=TRUE)
## png
## 2
plot_mri_style_image_human_all_slices(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/evoHumanHistory_regionwise_peak",transistionsToUseForFunctionalMapPostfix,".png"),evolutionaryAgeUniqueRegionsHumanAtlasWiseTop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=TRUE,isAlreadyInUniqueRegionsFormat=TRUE)
## png
## 2
plot_mri_style_image_human_all_slices(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/evoHumanHistory_peak_modeRegion",transistionsToUseForFunctionalMapPostfix,".png"),evolutionaryAgeBiopsyWiseTop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=TRUE,isAlreadyInUniqueRegionsFormat=TRUE)
## png
## 2
plot_mri_style_image_human_all_slices(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/evoHumanHistory_peak_modeRegion",transistionsToUseForFunctionalMapPostfix,"_sign.png"),evolutionaryAgeBiopsyWiseTopSign,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=TRUE,isAlreadyInUniqueRegionsFormat=TRUE)
## png
## 2
plot_mri_style_image_human_all_slices(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/evoHumanHistoryWA_meanRegion",transistionsToUseForFunctionalMapPostfix,".png"),evolutionaryAgeBiopsyWiseWATop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=FALSE,isAlreadyInUniqueRegionsFormat=TRUE)
## png
## 2
getColorScale<-function(amountOfColors){
return(colorpanel(amountOfColors,"#00007F","#F4F4F4","#7F0000"))
}
plot_mri_style_image_human(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistoryWA_meanRegion_br",transistionsToUseForFunctionalMapPostfix,".png"),evolutionaryAgeBiopsyWiseWATop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=FALSE,isAlreadyInUniqueRegionsFormat=TRUE)
## png
## 2
plot_mri_style_image_human_all_slices(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted_full/evoHumanHistoryWA_meanRegion_br",transistionsToUseForFunctionalMapPostfix,".png"),evolutionaryAgeBiopsyWiseWATop,length(transistionsToUseForFunctionalMap),mip=FALSE,discreteColorScale=FALSE,isAlreadyInUniqueRegionsFormat=TRUE)
## png
## 2
getColorScale<-function(amountOfColors){
return(colorpanel(amountOfColors,"#ffffb2","#feb24c","#bd0026"))
}
# maxVal<-quantile(humanFunctionalMapsLogPvalueTop[humanFunctionalMapsLogPvalueTop<Inf],0.95,na.rm=TRUE)
# for(actTrans in 1:length(transistionsToUseForFunctionalMap)){
# plot_mri_style_image_human(paste0("../storage/",prefix,"evoHistoryHuman_dndsweighted/evoHumanHistory_t",transistionsToUseForFunctionalMap[actTrans],"_mean_of_significant_sites.png"),humanFunctionalMapsLogPvalueTop[actTrans,],maxVal,mip=TRUE)
# }
#Calculate biclustering
############ RUN UNTIL HERE, THEN COMPUTE BICLUSTERING ############
############ YOU SHOULD BETTER RUN THIS ON A MACHINE ############
############ WITH >10 CORES SINCE THIS CAN TAKE A WHILE ############
#you don't need to run this, results are included
#source("GABI_BICLUSTERING.R")
for(actBiclustering in c("ALL_DNDS_","C_mouse_to_human_","E2_hominids_speciation_with_ancestorsNOG_","C3_mouse_to_human_late_no_gorgo_")){
for(lowOrHigh in c("high","low")){
if(lowOrHigh=="high"){
biclustering_results<-readRDS(paste0("../storage/biclustering_results/",actBiclustering,"ALL_NETWORKS_0_9_0_75.rds"))
}else{
biclustering_results<-readRDS(paste0("../storage/biclustering_results/",actBiclustering,"ALL_NETWORKS_0_1_0_75.rds"))
}
if(actBiclustering=="ALL_DNDS_"){
time_to_use <- c(1:10,11, 12, 13, 14, 15, 18, 19, 16, 17, 20, 21, 10)
#remove unstable clusters
if(lowOrHigh=="high"){
}else{
biclustering_results<-biclustering_results[-c(21)]
}
}
if(actBiclustering=="C_mouse_to_human_"){
time_to_use <- c(11, 12, 13, 14, 15, 18, 19, 16, 17, 20, 21, 10)
#remove unstable clusters
if(lowOrHigh=="high"){
}else{
biclustering_results<-biclustering_results[-c(21)]
}
}
if(actBiclustering=="E2_hominids_speciation_with_ancestorsNOG_"){
time_to_use <- c(18, 15, 7, 19, 16, 17, 20, 21, 8, 9, 10)
#remove unstable clusters
if(lowOrHigh=="high"){
biclustering_results<-biclustering_results[-c(6,10)]
}else{
biclustering_results<-biclustering_results[-c(21)]
}
}
if(actBiclustering=="C3_mouse_to_human_late_no_gorgo_"){
time_to_use <- c(10, 16, 17, 20, 21)
#remove unstable clusters
if(lowOrHigh=="high"){
}else{
biclustering_results<-biclustering_results[-c(15)]
}
}
biclusteringSummary <- c()
clusterFeatures <- lapply(biclustering_results, function(x) {
dndsCols <- x$features[x$features <= ncol(normalized_dnds_per_gene)]
networkCols <- x$features[x$features > ncol(normalized_dnds_per_gene)]
c(
paste(colnames(normalized_dnds_per_gene)[dndsCols], collapse = " "),
paste(colnames(gene_wise_network_correlation)[networkCols - ncol(normalized_dnds_per_gene)], collapse =
" "),
length(x$samples),
paste(apply(as.matrix(normalized_dnds_per_gene[x$samples, dndsCols]), 2, function(x) {
round(mean(x, na.rm = TRUE), digits = 3)
}), collapse = " "),
paste(apply(as.matrix(normalized_dnds_per_gene[x$samples, dndsCols]), 2, function(x) {
round(sd(x, na.rm = TRUE), digits = 3)
}), collapse = " "),
paste(apply(as.matrix(normalized_dnds_per_gene[x$samples, dndsCols]), 2, function(x) {
round(quantile(x, na.rm = TRUE, probs = 0.1), digits = 3)
}), collapse = " "),
paste(apply(as.matrix(normalized_dnds_per_gene[x$samples, dndsCols]), 2, function(x) {
round(quantile(x, na.rm = TRUE, probs = 0.9), digits = 3)
}), collapse = " "),
paste(apply(as.matrix(dnds_per_gene[x$samples, dndsCols]), 2, function(x) {
round(mean(x, na.rm = TRUE), digits = 3)
}), collapse = " "),
paste(apply(as.matrix(dnds_per_gene[x$samples, dndsCols]), 2, function(x) {
round(sd(x, na.rm = TRUE), digits = 3)
}), collapse = " "),
paste(apply(as.matrix(dnds_per_gene[x$samples, dndsCols]), 2, function(x) {
round(quantile(x, na.rm = TRUE, probs = 0.1), digits = 3)
}), collapse = " "),
paste(apply(as.matrix(dnds_per_gene[x$samples, dndsCols]), 2, function(x) {
round(quantile(x, na.rm = TRUE, probs = 0.9), digits = 3)
}), collapse = " ")
)
})
for (actCluster in 1:length(clusterFeatures)) {
biclusteringSummary <- rbind(biclusteringSummary, clusterFeatures[[actCluster]])
}
colnames(biclusteringSummary) <-
c(
"DNDS columns",
"Network columns",
"Amount of Genes",
"RAW DNDS",
"STD DNDS",
"0.1 Quantile DNDS",
"0.9 Quantile DNDS",
"Mean RAW DNDS",
"STD RAW DNDS",
"0.1 Quantile RAW DNDS",
"0.9 RAW Quantile DNDS"
)
rownames(biclusteringSummary) <-
paste0("Bicluster_", c(1:nrow(biclusteringSummary)))
print(kable(biclusteringSummary,caption=gsub("_", " ", paste0(actBiclustering,"_",lowOrHigh))))
}
}
ALL DNDS high
| Bicluster_1 |
HUMAN |
LA-AVG |
21 |
0.96 |
0.03 |
0.916 |
0.995 |
1.143 |
0.432 |
0.77 |
1.719 |
| Bicluster_2 |
HUMAN |
WM-REST |
19 |
0.962 |
0.03 |
0.915 |
0.995 |
1.178 |
0.44 |
0.768 |
1.745 |
| Bicluster_3 |
HUMAN |
MO-RH+LH |
18 |
0.959 |
0.026 |
0.926 |
0.989 |
1.065 |
0.377 |
0.789 |
1.478 |
| Bicluster_4 |
HUMAN |
WM-REST WM-FACE |
14 |
0.963 |
0.033 |
0.914 |
0.996 |
1.235 |
0.483 |
0.765 |
1.809 |
| Bicluster_5 |
HUMAN |
cortico-limbic_network |
13 |
0.954 |
0.028 |
0.926 |
0.992 |
1.04 |
0.366 |
0.787 |
1.62 |
| Bicluster_6 |
MOCMNGPH_ND |
deafult_mode_network |
8 |
0.949 |
0.036 |
0.912 |
0.988 |
0.945 |
0.196 |
0.758 |
1.176 |
| Bicluster_7 |
MOUSE OTOGA |
salience_network |
93 |
0.957 0.955 |
0.027 0.028 |
0.92 0.915 |
0.993 0.995 |
0.42 0.559 |
0.154 0.2 |
0.302 0.388 |
0.57 0.881 |
| Bicluster_8 |
NEANDERTHAL |
GA-AVG |
7 |
0.93 |
0.023 |
0.903 |
0.954 |
0.815 |
0.085 |
0.719 |
0.915 |
| Bicluster_9 |
HUMAN |
GA-REWARD WM-REST |
11 |
0.961 |
0.031 |
0.916 |
0.995 |
1.161 |
0.454 |
0.77 |
1.719 |
| Bicluster_10 |
HUMAN |
EM-FACES MO-RH+LH |
9 |
0.963 |
0.031 |
0.915 |
0.988 |
1.142 |
0.337 |
0.768 |
1.452 |
| Bicluster_11 |
MOCMNGPH_ND |
sensorimotor_network |
5 |
0.945 |
0.039 |
0.91 |
0.986 |
0.923 |
0.209 |
0.751 |
1.154 |
| Bicluster_12 |
MOCMNDHN_PG |
visuospatial_network |
4 |
0.958 |
0.03 |
0.935 |
0.987 |
1.274 |
0.556 |
0.915 |
1.797 |
| Bicluster_13 |
OTOGA |
mPFC-AMY_network mPFC-Acb_network |
31 |
0.944 |
0.031 |
0.91 |
0.989 |
0.507 |
0.174 |
0.377 |
0.718 |
| Bicluster_14 |
MOCMNGPH_ND |
EM-FACES MO-RH+LH RE-AVG |
4 |
0.948 |
0.035 |
0.923 |
0.982 |
1.023 |
0.419 |
0.796 |
1.408 |
| Bicluster_15 |
DENISOVAN |
fronto-parietal |
44 |
0.957 |
0.03 |
0.918 |
0.994 |
1.884 |
1.042 |
1.172 |
2.937 |
| Bicluster_16 |
MOCMNGP_HND |
central_executive_network |
27 |
0.948 |
0.034 |
0.907 |
0.996 |
1.393 |
0.711 |
0.834 |
2.549 |
| Bicluster_17 |
DENISOVAN |
ventral_attention |
37 |
0.954 |
0.026 |
0.923 |
0.99 |
1.663 |
0.568 |
1.198 |
2.479 |
| Bicluster_18 |
NOMLE |
MO-RF+LF |
26 |
0.947 |
0.029 |
0.906 |
0.98 |
0.869 |
0.266 |
0.597 |
1.15 |
| Bicluster_19 |
DENISOVAN |
dorsal_attention |
16 |
0.95 |
0.027 |
0.922 |
0.983 |
1.594 |
0.543 |
1.194 |
2.099 |
| Bicluster_20 |
PANTR |
SO-TOM |
7 |
0.958 |
0.036 |
0.914 |
0.993 |
1.446 |
0.7 |
0.85 |
2.154 |
ALL DNDS low
| Bicluster_1 |
MOCMNGPD_HN |
LA-AVG |
281 |
0.001 |
0.006 |
0 |
0 |
0.002 |
0.013 |
0 |
0 |
| Bicluster_2 |
MOCMNGPD_HN |
WM-REST |
222 |
0.001 |
0.008 |
0 |
0 |
0.002 |
0.016 |
0 |
0 |
| Bicluster_3 |
MOCMNGPD_HN |
cortico-limbic_network |
215 |
0.001 |
0.007 |
0 |
0 |
0.002 |
0.015 |
0 |
0 |
| Bicluster_4 |
MOCMNGPD_HN |
deafult_mode_network |
199 |
0 |
0.005 |
0 |
0 |
0.001 |
0.01 |
0 |
0 |
| Bicluster_5 |
MOCMNGPD_HN |
salience_network |
179 |
0.003 |
0.014 |
0 |
0 |
0.005 |
0.025 |
0 |
0 |
| Bicluster_6 |
MOCMNGPH_ND |
WM-FACE |
179 |
0.001 |
0.008 |
0 |
0 |
0.002 |
0.017 |
0 |
0 |
| Bicluster_7 |
MOCMNGPH_ND |
GA-AVG |
181 |
0.001 |
0.008 |
0 |
0 |
0.003 |
0.018 |
0 |
0 |
| Bicluster_8 |
MOCMNGPH_ND |
MO-RH+LH |
160 |
0.001 |
0.01 |
0 |
0 |
0.003 |
0.018 |
0 |
0 |
| Bicluster_9 |
MOCMNGPH_ND |
EM-FACES |
105 |
0.001 |
0.01 |
0 |
0 |
0.002 |
0.015 |
0 |
0 |
| Bicluster_10 |
MOCMNGPD_HN |
GA-REWARD |
73 |
0.002 |
0.012 |
0 |
0 |
0.005 |
0.024 |
0 |
0 |
| Bicluster_11 |
MOCMNGPD_HN |
mPFC-AMY_network |
66 |
0.002 |
0.011 |
0 |
0 |
0.004 |
0.021 |
0 |
0 |
| Bicluster_12 |
MOCMNGPD_HN |
visuospatial_network |
64 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_13 |
MOCMNGPD_HN |
sensorimotor_network |
64 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_14 |
MOCMNGPD_HN |
fronto-parietal |
62 |
0 |
0.004 |
0 |
0 |
0.002 |
0.013 |
0 |
0 |
| Bicluster_15 |
MOCMNGPH_ND |
central_executive_network |
70 |
0.002 |
0.011 |
0 |
0 |
0.005 |
0.023 |
0 |
0 |
| Bicluster_16 |
MOCMNGPD_HN |
mPFC-Acb_network |
54 |
0.001 |
0.008 |
0 |
0 |
0.002 |
0.016 |
0 |
0 |
| Bicluster_17 |
MOCMNGPH_ND |
RE-AVG |
51 |
0.002 |
0.013 |
0 |
0 |
0.003 |
0.02 |
0 |
0 |
| Bicluster_18 |
MOCMNGPD_HN |
ventral_attention |
37 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_19 |
NEANDERTHAL |
MO-RF+LF |
31 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_20 |
MOCMNGPD_HN |
dorsal_attention |
19 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_21 |
MOUSE OTOGA GORGO MO_CMNGPDNH MOC_MNGPDNH MOCMN_GPDNH |
LA-AVG MO-RH+LH MO-T |
4 |
0.03 0.009 0 0 0 0 |
0.034 0.017 0 0 0 0.001 |
0.01 0 0 0 0 0 |
0.061 0.024 0 0 0 0.001 |
0.009 0.003 0 0 0 0 |
0.007 0.007 0 0 0 0 |
0.005 0 0 0 0 0 |
0.016 0.01 0 0 0 0 |
C mouse to human high
| Bicluster_1 |
HUMAN |
LA-AVG |
21 |
0.96 |
0.03 |
0.916 |
0.995 |
1.143 |
0.432 |
0.77 |
1.719 |
| Bicluster_2 |
HUMAN |
WM-REST |
19 |
0.962 |
0.03 |
0.915 |
0.995 |
1.178 |
0.44 |
0.768 |
1.745 |
| Bicluster_3 |
HUMAN |
MO-RH+LH |
18 |
0.959 |
0.026 |
0.926 |
0.989 |
1.065 |
0.377 |
0.789 |
1.478 |
| Bicluster_4 |
HUMAN |
WM-REST WM-FACE |
14 |
0.963 |
0.033 |
0.914 |
0.996 |
1.235 |
0.483 |
0.765 |
1.809 |
| Bicluster_5 |
MOCMNGPH_ND |
cortico-limbic_network |
8 |
0.955 |
0.033 |
0.912 |
0.988 |
0.967 |
0.182 |
0.758 |
1.176 |
| Bicluster_6 |
MOCMNGPH_ND |
deafult_mode_network |
8 |
0.949 |
0.036 |
0.912 |
0.988 |
0.945 |
0.196 |
0.758 |
1.176 |
| Bicluster_7 |
MOCMN_GPDNH |
salience_network |
129 |
0.949 |
0.028 |
0.91 |
0.985 |
1.25 |
1.172 |
0.713 |
1.653 |
| Bicluster_8 |
HUMAN |
GA-AVG LA-AVG |
11 |
0.959 |
0.026 |
0.94 |
0.992 |
1.065 |
0.356 |
0.817 |
1.619 |
| Bicluster_9 |
HUMAN |
GA-REWARD WM-REST |
11 |
0.961 |
0.031 |
0.916 |
0.995 |
1.161 |
0.454 |
0.77 |
1.719 |
| Bicluster_10 |
HUMAN |
EM-FACES MO-RH+LH |
9 |
0.963 |
0.031 |
0.915 |
0.988 |
1.142 |
0.337 |
0.768 |
1.452 |
| Bicluster_11 |
MOCMNGPH_ND |
sensorimotor_network |
5 |
0.945 |
0.039 |
0.91 |
0.986 |
0.923 |
0.209 |
0.751 |
1.154 |
| Bicluster_12 |
MOCMNDHN_PG |
visuospatial_network |
4 |
0.958 |
0.03 |
0.935 |
0.987 |
1.274 |
0.556 |
0.915 |
1.797 |
| Bicluster_13 |
MOCMNGPH_ND |
EM-FACES MO-RH+LH RE-AVG |
4 |
0.948 |
0.035 |
0.923 |
0.982 |
1.023 |
0.419 |
0.796 |
1.408 |
| Bicluster_14 |
MOCMNGP_HND |
central_executive_network |
27 |
0.948 |
0.034 |
0.907 |
0.996 |
1.393 |
0.711 |
0.834 |
2.549 |
| Bicluster_15 |
MOC_MNGPDNH |
mPFC-Acb_network |
37 |
0.943 |
0.03 |
0.905 |
0.983 |
1.109 |
0.692 |
0.668 |
1.622 |
| Bicluster_16 |
MOC_MNGPDNH |
mPFC-AMY_network |
33 |
0.948 |
0.032 |
0.91 |
0.993 |
1.284 |
0.902 |
0.691 |
2.356 |
| Bicluster_17 |
MOC_MNGPDNH |
fronto-parietal |
30 |
0.948 |
0.029 |
0.906 |
0.987 |
1.154 |
0.698 |
0.673 |
1.799 |
| Bicluster_18 |
MOCMNGPH_ND |
LA-AVG cortico-limbic_network ventral_attention |
3 |
0.954 |
0.019 |
0.938 |
0.966 |
0.915 |
0.094 |
0.837 |
0.983 |
| Bicluster_19 |
MO_CMNGPDNH |
MO-RF+LF |
24 |
0.953 |
0.029 |
0.912 |
0.987 |
0.694 |
0.308 |
0.451 |
0.895 |
| Bicluster_20 |
MOCM_NGPDNH |
dorsal_attention |
9 |
0.967 |
0.021 |
0.946 |
0.99 |
1.503 |
0.529 |
1.03 |
2.34 |
| Bicluster_21 |
MOCM_NGPDNH |
MO-RH+LH SO-TOM WM-REST |
5 |
0.952 |
0.04 |
0.912 |
0.992 |
2.698 |
3.594 |
0.78 |
6.144 |
C mouse to human low
| Bicluster_1 |
MOCMNGPD_HN |
LA-AVG |
281 |
0.001 |
0.006 |
0 |
0 |
0.002 |
0.013 |
0 |
0 |
| Bicluster_2 |
MOCMNGPD_HN |
WM-REST |
222 |
0.001 |
0.008 |
0 |
0 |
0.002 |
0.016 |
0 |
0 |
| Bicluster_3 |
MOCMNGPD_HN |
cortico-limbic_network |
215 |
0.001 |
0.007 |
0 |
0 |
0.002 |
0.015 |
0 |
0 |
| Bicluster_4 |
MOCMNGPD_HN |
deafult_mode_network |
199 |
0 |
0.005 |
0 |
0 |
0.001 |
0.01 |
0 |
0 |
| Bicluster_5 |
MOCMNGPD_HN |
salience_network |
179 |
0.003 |
0.014 |
0 |
0 |
0.005 |
0.025 |
0 |
0 |
| Bicluster_6 |
MOCMNGPH_ND |
WM-FACE |
179 |
0.001 |
0.008 |
0 |
0 |
0.002 |
0.017 |
0 |
0 |
| Bicluster_7 |
MOCMNGPH_ND |
GA-AVG |
181 |
0.001 |
0.008 |
0 |
0 |
0.003 |
0.018 |
0 |
0 |
| Bicluster_8 |
MOCMNGPH_ND |
MO-RH+LH |
160 |
0.001 |
0.01 |
0 |
0 |
0.003 |
0.018 |
0 |
0 |
| Bicluster_9 |
MOCMNGPH_ND |
EM-FACES |
105 |
0.001 |
0.01 |
0 |
0 |
0.002 |
0.015 |
0 |
0 |
| Bicluster_10 |
MOCMNGPD_HN |
GA-REWARD |
73 |
0.002 |
0.012 |
0 |
0 |
0.005 |
0.024 |
0 |
0 |
| Bicluster_11 |
MOCMNGPD_HN |
mPFC-AMY_network |
66 |
0.002 |
0.011 |
0 |
0 |
0.004 |
0.021 |
0 |
0 |
| Bicluster_12 |
MOCMNGPH_ND |
visuospatial_network |
86 |
0.002 |
0.01 |
0 |
0 |
0.004 |
0.021 |
0 |
0 |
| Bicluster_13 |
MOCMNGPD_HN |
sensorimotor_network |
64 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_14 |
MOCMNGPD_HN |
fronto-parietal |
62 |
0 |
0.004 |
0 |
0 |
0.002 |
0.013 |
0 |
0 |
| Bicluster_15 |
MOCMNGPD_HN |
central_executive_network |
58 |
0.003 |
0.013 |
0 |
0 |
0.006 |
0.026 |
0 |
0 |
| Bicluster_16 |
MOCMNGPD_HN |
mPFC-Acb_network |
54 |
0.001 |
0.008 |
0 |
0 |
0.002 |
0.016 |
0 |
0 |
| Bicluster_17 |
MOCMNGPH_ND |
RE-AVG |
51 |
0.002 |
0.013 |
0 |
0 |
0.003 |
0.02 |
0 |
0 |
| Bicluster_18 |
MOCMNGPD_HN |
ventral_attention |
37 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_19 |
MOCMNGPH_ND |
MO-RF+LF |
34 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_20 |
MOCMNGPD_HN |
dorsal_attention |
19 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_21 |
MO_CMNGPDNH MOC_MNGPDNH MOCM_NGPDNH MOCMN_GPDNH |
LA-AVG MO-RH+LH MO-T |
4 |
0 0 0 0 |
0 0 0 0.001 |
0 0 0 0 |
0 0 0 0.001 |
0 0 0 0 |
0 0 0 0 |
0 0 0 0 |
0 0 0 0 |
E2 hominids speciation with ancestorsNOG high
| Bicluster_1 |
HUMAN |
LA-AVG |
21 |
0.96 |
0.03 |
0.916 |
0.995 |
1.143 |
0.432 |
0.77 |
1.719 |
| Bicluster_2 |
NEANDERTHAL |
WM-REST |
10 |
0.939 |
0.026 |
0.907 |
0.962 |
0.87 |
0.145 |
0.735 |
0.957 |
| Bicluster_3 |
HUMAN |
MO-RH+LH |
18 |
0.959 |
0.026 |
0.926 |
0.989 |
1.065 |
0.377 |
0.789 |
1.478 |
| Bicluster_4 |
DENISOVAN |
salience_network |
188 |
0.946 |
0.03 |
0.907 |
0.99 |
1.634 |
0.742 |
1.122 |
2.518 |
| Bicluster_5 |
DENISOVAN |
cortico-limbic_network deafult_mode_network |
106 |
0.945 |
0.027 |
0.91 |
0.983 |
1.543 |
0.584 |
1.14 |
2.07 |
| Bicluster_6 |
NEANDERTHAL |
GA-AVG |
7 |
0.93 |
0.023 |
0.903 |
0.954 |
0.815 |
0.085 |
0.719 |
0.915 |
| Bicluster_7 |
HUMAN |
GA-REWARD WM-REST |
11 |
0.961 |
0.031 |
0.916 |
0.995 |
1.161 |
0.454 |
0.77 |
1.719 |
| Bicluster_8 |
HUMAN |
EM-FACES MO-RH+LH |
9 |
0.963 |
0.031 |
0.915 |
0.988 |
1.142 |
0.337 |
0.768 |
1.452 |
| Bicluster_9 |
DENISOVAN |
sensorimotor_network |
62 |
0.951 |
0.028 |
0.916 |
0.989 |
1.731 |
0.901 |
1.164 |
2.442 |
| Bicluster_10 |
DENISOVAN |
mPFC-Acb_network |
45 |
0.951 |
0.026 |
0.919 |
0.984 |
1.628 |
0.742 |
1.176 |
2.123 |
| Bicluster_11 |
MOCMNGPH_ND |
EM-FACES MO-RH+LH RE-AVG |
4 |
0.948 |
0.035 |
0.923 |
0.982 |
1.023 |
0.419 |
0.796 |
1.408 |
| Bicluster_12 |
DENISOVAN |
fronto-parietal |
44 |
0.957 |
0.03 |
0.918 |
0.994 |
1.884 |
1.042 |
1.172 |
2.937 |
| Bicluster_13 |
MOCMNGP_HND |
central_executive_network |
27 |
0.948 |
0.034 |
0.907 |
0.996 |
1.393 |
0.711 |
0.834 |
2.549 |
| Bicluster_14 |
DENISOVAN |
mPFC-AMY_network |
39 |
0.944 |
0.028 |
0.908 |
0.981 |
1.561 |
0.788 |
1.126 |
2.029 |
| Bicluster_15 |
DENISOVAN |
ventral_attention |
37 |
0.954 |
0.026 |
0.923 |
0.99 |
1.663 |
0.568 |
1.198 |
2.479 |
| Bicluster_16 |
DENISOVAN |
MO-RF+LF MO-RH+LH |
27 |
0.953 |
0.028 |
0.915 |
0.992 |
1.672 |
0.586 |
1.157 |
2.641 |
| Bicluster_17 |
DENISOVAN |
dorsal_attention |
16 |
0.95 |
0.027 |
0.922 |
0.983 |
1.594 |
0.543 |
1.194 |
2.099 |
| Bicluster_18 |
PANTR |
SO-TOM |
7 |
0.958 |
0.036 |
0.914 |
0.993 |
1.446 |
0.7 |
0.85 |
2.154 |
E2 hominids speciation with ancestorsNOG low
| Bicluster_1 |
MOCMNGPD_HN |
LA-AVG |
281 |
0.001 |
0.006 |
0 |
0 |
0.002 |
0.013 |
0 |
0 |
| Bicluster_2 |
MOCMNGPD_HN |
WM-REST |
222 |
0.001 |
0.008 |
0 |
0 |
0.002 |
0.016 |
0 |
0 |
| Bicluster_3 |
MOCMNGPD_HN |
cortico-limbic_network |
215 |
0.001 |
0.007 |
0 |
0 |
0.002 |
0.015 |
0 |
0 |
| Bicluster_4 |
NEANDERTHAL |
deafult_mode_network |
190 |
0.004 |
0.017 |
0 |
0 |
0.005 |
0.025 |
0 |
0 |
| Bicluster_5 |
MOCMNGPD_HN |
salience_network |
179 |
0.003 |
0.014 |
0 |
0 |
0.005 |
0.025 |
0 |
0 |
| Bicluster_6 |
MOCMNGPD_HN |
WM-FACE |
134 |
0.001 |
0.006 |
0 |
0 |
0.002 |
0.014 |
0 |
0 |
| Bicluster_7 |
NEANDERTHAL |
GA-AVG |
122 |
0.003 |
0.014 |
0 |
0 |
0.006 |
0.024 |
0 |
0 |
| Bicluster_8 |
NEANDERTHAL |
MO-RH+LH |
115 |
0.002 |
0.011 |
0 |
0 |
0.003 |
0.018 |
0 |
0 |
| Bicluster_9 |
NEANDERTHAL |
EM-FACES |
81 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_10 |
MOCMNGPD_HN |
GA-REWARD |
73 |
0.002 |
0.012 |
0 |
0 |
0.005 |
0.024 |
0 |
0 |
| Bicluster_11 |
MOCMNGPD_HN |
mPFC-AMY_network |
66 |
0.002 |
0.011 |
0 |
0 |
0.004 |
0.021 |
0 |
0 |
| Bicluster_12 |
MOCMNGPD_HN |
visuospatial_network |
64 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_13 |
MOCMNGPD_HN |
sensorimotor_network |
64 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_14 |
NEANDERTHAL |
central_executive_network |
56 |
0.001 |
0.007 |
0 |
0 |
0.003 |
0.015 |
0 |
0 |
| Bicluster_15 |
MOCMNGPD_HN |
fronto-parietal |
62 |
0 |
0.004 |
0 |
0 |
0.002 |
0.013 |
0 |
0 |
| Bicluster_16 |
MOCMNGPD_HN |
mPFC-Acb_network |
54 |
0.001 |
0.008 |
0 |
0 |
0.002 |
0.016 |
0 |
0 |
| Bicluster_17 |
MOCMNGPD_HN |
RE-AVG |
40 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_18 |
MOCMNGPD_HN |
ventral_attention |
37 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_19 |
NEANDERTHAL |
MO-RF+LF |
31 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_20 |
MOCMNGPD_HN |
dorsal_attention |
19 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
C3 mouse to human late no gorgo high
| Bicluster_1 |
MOCMNGP_HND |
LA-AVG |
102 |
0.944 |
0.029 |
0.908 |
0.987 |
1.247 |
0.551 |
0.838 |
1.893 |
| Bicluster_2 |
MOCMNGP_HND |
cortico-limbic_network |
85 |
0.948 |
0.031 |
0.909 |
0.995 |
1.354 |
0.647 |
0.846 |
2.353 |
| Bicluster_3 |
MOCMNGP_HND |
WM-REST |
79 |
0.947 |
0.029 |
0.908 |
0.985 |
1.283 |
0.554 |
0.837 |
1.767 |
| Bicluster_4 |
HUMAN |
MO-RH+LH |
18 |
0.959 |
0.026 |
0.926 |
0.989 |
1.065 |
0.377 |
0.789 |
1.478 |
| Bicluster_5 |
MOCMNGP_HND |
deafult_mode_network |
72 |
0.948 |
0.031 |
0.91 |
0.992 |
1.327 |
0.612 |
0.854 |
2.117 |
| Bicluster_6 |
MOCMNGP_HND |
salience_network |
71 |
0.953 |
0.026 |
0.921 |
0.989 |
1.313 |
0.44 |
0.915 |
1.982 |
| Bicluster_7 |
MOCMNGP_HND |
WM-FACE |
54 |
0.951 |
0.027 |
0.914 |
0.984 |
1.326 |
0.579 |
0.876 |
1.715 |
| Bicluster_8 |
HUMAN |
GA-REWARD WM-REST |
11 |
0.961 |
0.031 |
0.916 |
0.995 |
1.161 |
0.454 |
0.77 |
1.719 |
| Bicluster_9 |
HUMAN |
GA-AVG LA-AVG |
11 |
0.959 |
0.026 |
0.94 |
0.992 |
1.065 |
0.356 |
0.817 |
1.619 |
| Bicluster_10 |
MOCMNGP_HND |
EM-FACES |
34 |
0.956 |
0.027 |
0.917 |
0.984 |
1.349 |
0.481 |
0.896 |
1.728 |
| Bicluster_11 |
MOCMNGP_HND |
central_executive_network |
27 |
0.948 |
0.034 |
0.907 |
0.996 |
1.393 |
0.711 |
0.834 |
2.549 |
| Bicluster_12 |
MOCMNGP_HND |
visuospatial_network |
27 |
0.948 |
0.032 |
0.904 |
0.986 |
1.265 |
0.414 |
0.821 |
1.833 |
| Bicluster_13 |
MOCMNGP_HND |
sensorimotor_network |
24 |
0.943 |
0.029 |
0.904 |
0.985 |
1.215 |
0.51 |
0.824 |
1.832 |
| Bicluster_14 |
MOCMNGP_HND |
mPFC-AMY_network |
22 |
0.945 |
0.028 |
0.913 |
0.988 |
1.297 |
0.71 |
0.874 |
1.92 |
| Bicluster_15 |
MOCMNGP_HND |
mPFC-Acb_network |
20 |
0.952 |
0.028 |
0.916 |
0.995 |
1.407 |
0.743 |
0.887 |
2.42 |
| Bicluster_16 |
MOCMNGP_HND |
fronto-parietal |
19 |
0.948 |
0.03 |
0.907 |
0.981 |
1.295 |
0.546 |
0.837 |
1.736 |
| Bicluster_17 |
MOCMNGP_HND |
RE-AVG |
19 |
0.956 |
0.033 |
0.913 |
0.994 |
1.484 |
0.685 |
0.868 |
2.299 |
| Bicluster_18 |
MOCMNGP_HND |
MO-RF+LF |
16 |
0.962 |
0.029 |
0.926 |
0.994 |
1.578 |
0.779 |
0.942 |
2.662 |
| Bicluster_19 |
MOCMNGPH_ND |
LA-AVG cortico-limbic_network ventral_attention |
3 |
0.954 |
0.019 |
0.938 |
0.966 |
0.915 |
0.094 |
0.837 |
0.983 |
| Bicluster_20 |
MOCMNGP_HND |
SO-TOM WM-REST |
3 |
0.933 |
0.03 |
0.911 |
0.958 |
1.033 |
0.256 |
0.858 |
1.245 |
| Bicluster_21 |
MOCMNGP_HND |
fronto-parietal dorsal_attention |
3 |
0.953 |
0.031 |
0.927 |
0.973 |
1.248 |
0.318 |
0.983 |
1.478 |
C3 mouse to human late no gorgo low
| Bicluster_1 |
HUMAN |
LA-AVG |
493 |
0.001 |
0.008 |
0 |
0 |
0.003 |
0.016 |
0 |
0 |
| Bicluster_2 |
HUMAN |
WM-REST |
395 |
0.001 |
0.006 |
0 |
0 |
0.002 |
0.014 |
0 |
0 |
| Bicluster_3 |
HUMAN |
cortico-limbic_network |
357 |
0.002 |
0.01 |
0 |
0 |
0.003 |
0.017 |
0 |
0 |
| Bicluster_4 |
HUMAN |
deafult_mode_network |
341 |
0.001 |
0.01 |
0 |
0 |
0.003 |
0.017 |
0 |
0 |
| Bicluster_5 |
HUMAN |
salience_network |
292 |
0.002 |
0.01 |
0 |
0 |
0.003 |
0.018 |
0 |
0 |
| Bicluster_6 |
HUMAN |
WM-FACE |
248 |
0 |
0.004 |
0 |
0 |
0.001 |
0.011 |
0 |
0 |
| Bicluster_7 |
HUMAN |
GA-AVG |
242 |
0.001 |
0.007 |
0 |
0 |
0.003 |
0.017 |
0 |
0 |
| Bicluster_8 |
HUMAN |
MO-RH+LH |
229 |
0.001 |
0.005 |
0 |
0 |
0.002 |
0.012 |
0 |
0 |
| Bicluster_9 |
HUMAN |
EM-FACES |
144 |
0.001 |
0.005 |
0 |
0 |
0.002 |
0.013 |
0 |
0 |
| Bicluster_10 |
HUMAN |
GA-REWARD |
129 |
0.002 |
0.008 |
0 |
0 |
0.004 |
0.019 |
0 |
0 |
| Bicluster_11 |
HUMAN |
visuospatial_network |
114 |
0.002 |
0.012 |
0 |
0 |
0.005 |
0.022 |
0 |
0 |
| Bicluster_12 |
HUMAN |
sensorimotor_network |
110 |
0.002 |
0.012 |
0 |
0 |
0.005 |
0.022 |
0 |
0 |
| Bicluster_13 |
HUMAN |
central_executive_network |
99 |
0.001 |
0.007 |
0 |
0 |
0.001 |
0.012 |
0 |
0 |
| Bicluster_14 |
HUMAN |
mPFC-AMY_network |
91 |
0.002 |
0.014 |
0 |
0 |
0.003 |
0.019 |
0 |
0 |
| Bicluster_15 |
HUMAN |
mPFC-Acb_network |
82 |
0.002 |
0.014 |
0 |
0 |
0.003 |
0.021 |
0 |
0 |
| Bicluster_16 |
HUMAN |
RE-AVG |
73 |
0.001 |
0.004 |
0 |
0 |
0.002 |
0.014 |
0 |
0 |
| Bicluster_17 |
HUMAN |
ventral_attention |
60 |
0.001 |
0.009 |
0 |
0 |
0.002 |
0.015 |
0 |
0 |
| Bicluster_18 |
HUMAN |
MO-RF+LF |
46 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_19 |
HUMAN |
fronto-parietal dorsal_attention |
28 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Bicluster_20 |
MOCMNGP_HND |
SO-TOM |
37 |
0.007 |
0.017 |
0 |
0.029 |
0.009 |
0.021 |
0 |
0.048 |